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ABSTRACT 

This  thesis  contains  discussion,  theory  and  program  code 
for  a  computational  fluid  dynamics  (CFD)  model  of  a  wing  of 
arbitrary  planform.  The  model  assumes  incompressible, 
inviscid,  irrotational  flow.  The  program  computes  forces 
acting  on  the  wing  by  modeling  the  flow  with  a  set  of  horse 
shoe  vortex  elements.  It  models  the  flow  over  an  arbitrary 
wing  using  two  solutions.  One  solution  is  the  ideal  lift, 
associated  with  a  cambered  and  twisted  wing.  The  other 
solution  is  the  additional  lift  associated  with  a  flat  wing. 
The  program  computes  wing  camber  and  twist  using  an  elliptic 
loading  distribution.  The  thesis  includes  the  FORTRAN  source 
code,  a  separable  User's  Manual  for  the  VORTEX  program, 
discussion  of  the  theory  applied  in  the  model,  and 
instructions  for  operating  the  program.  It  shows  a  sample 
wing  planform  with  tabular  and  graphic  results. 

The  thesis  also  discusses  two  other  CFD  models  based  on 
circulation  (T)  and  pressure  difference  (ACp).  It  presents 
some  of  the  problems  and  solutions  in  grid  generation. 
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I .  INTRODUCTION.  PURPOSE  AND  GOALS 

AE2035  is  an  introductory  course  in  aerodynamics  for 
aeronautical  engineering  students  at  Naval  Postgraduate 
School,  Monterey,  CA .  The  course  must  provide  students  with 
the  required  knowledge  for  future  courses  in  the 
Aeronautical  Engineering  Curriculum.  The  student  must 
understand  the  techniques  used  in  computational  fluid 
dynamics  (CFD)  to  satisfy  that  requirement.  Powerful 
computational  models  are  available  for  the  working 
ae r o dy nam i c i s t ,  but  they  are  seldom  usable  in  an 
introductory  course.  The  programming  language  for  these 
models  is  usually  optimized  for  computational  efficiency 
rather  than  teaching,  and  the  documentation  may  be  poor  or 
unavailable.  Students  need  a  computer  program  with  graphic 
output  which  is  less  complicated  to  use  as  a  teaching  tool. 
Providing  that  computer  program  is  the  goal  of  the  thesis. 

The  program  runs  on  a  micro  or  desk  top  computer,  in 
keeping  with  the  emphasis  on  use  of  individual  workstations 
to  supplement  the  school's  mainframe  computer.  The  micro 
computer  has  better  graphic  output  programs  available, 
without  the  extensive  programming  required  on  the  school 
mainframe.  This  flexibility  allows  the  student  to  modify  the 
output  to  suit  his/her  requirements. 


The  thesis  contains  five  principal  sections.  The  first 
section,  TECHNIQUES,  provides  a  brief  synopsis  of  three 
different  flow  models.  The  section  also  provides  some 
insight  to  certain  problems  that  arise  in  CFD. 

The  second  section,  VORTEX  PROGRAM  DEVELOPMENT,  contains 
a  detailed  discussion  of  the  horse  shoe  vortex  model.  The 
section  includes  the  theory  and  techniques  used  in  the 
VORTEX  program. 

The  third  section,  CIRCULATION  MODEL,  contains 
discussion  of  a  circulation  (T)  model.  The  section  includes 
the  theory  behind  the  model  and  some  of  the  problems 
enc  ount e  red. 

The  fourth  section,  PRESSURE  DIFFERENCE  MODEL,  contains 
discussion  of  the  pressure  difference  (ACp)  model.  The 
section  also  includes  the  theory  behind  the  model  and 
problems  encountered. 

The  fifth  and  last  principal  section,  the  USERS  MANUAL, 
is  separable  from  the  body  of  the  thesis  and  provides 
detailed  theory  on  the  horse  shoe  vortex  model.  The  section 
also  includes  instruction  on  use  of  the  program  and  an 
example  wing  planform  with  the  resulting  data. 

The  computer  listing  for  the  VORTEX  program  is  included 
in  the  Appendix. 


A.    PROGRAM  ASSUMPTIONS  AND  REQUIREMENTS 

To  keep  the  model  streamlined,  the  following  assumptions 
were  made  and  requirements  established. 

1.  Flow  is  steady,  inviscid,  i r r o t a t i ona  1  ,  and 
incompressible . 

2.  The  wing  surface  has  zero  thickness. 

3.  The  planform  has  any  sweep,  taper,  and  aspect  ratio. 

4.  The  program  must  support  low  aspect  ratio  planforms 
(less  than  1.0). 

5.  The  program  must  contain  comments,  for  easy 
modification  and  change. 

6.  The  program  results  must  be  available  in  tabular  and 
graphic  form . 

7.  The  program  must  use  accepted  variable  names,  (aspect 
ratio,  taper  ratio,  etc.) 

8.  The  program  must  have  two  options.  One  is  to  generate 
loading  for  a  flat  wing  shape.  The  other  is  to 
generate  a  wing  shape  for  an  elliptic  loading. 

9.  The  programming  language  must  be  FORTRAN,  with 
executable  and  source  files  provided. 


II  .  TECHNIQUES 

Three  different  models  were  investigated.  Each  of  the 
models  uses  a  governing  equation  for  the  induced  velocity  at 
each  control  point-1-  on  the  wing.  That  velocity  is 
associated  with  the  vortex  strengths  located  at  field  points 
on  the  wing.  The  induced  velocity  is  a  function  of  two 
factors.  One  is  the  vortex  strength  and  the  other  is  the 
physical  orientation  of  the  field  and  control  points.  Flow 
must  be  tangent  to  the  wing  surface  and  is  found  by  adding 
the  induced  velocity  vector  (w)  to  the  remote  velocity 
vector  (Vco).  The  end  result  is  a  direct  correlation  between 
the  wing's  shape'1  and  the  strengths  of  the  vortices.  This 
relationship  allows  either  the  strengths  of  the  vortices  or 
the  wing  shape  to  be  the  independent  variables. 


-'-Field  points  are  the  positions  of  the  vortex  elements, 
and  control  points  are  the  positions  where  the  velocity  is 
evaluated.  For  example,  a  control  point  is  where  flow 
tangency  is  enforced. 

^Wing  shape  will  be  used  throughout  to  describe  the 
surface  slope  of  the  wing.  For  example,  one  wing  shape 
might  be  a  flat  plate,  another  might  have  the  tip  twisted 
relative  to  the  root  chord. 


The  three  models  use  circulation  (T),  incremental 
circulation  (AT),  and  pressure  difference  coefficient  (ACp). 
These  quantities  are  related  as  follows. 

V  Ax  J  \dx  I  2  1 

AC  =  =  z  •  L 

p  V  V 

A.   HORSE  SHOE  VORTEX  MODEL 

This  program  models  the  flow  by  a  set  of  horse  shoe 
vortices  distributed  over  the  wing.  Each  vortex  consists  of 
a  bound  portion  perpendicular  to  the  remote  velocity,  plus 
two  trailing  portions.  Each  vortex  element  generates  an 
incremental  circulation  (AT)  around  the  element. 

The  program  generates  a  matrix  of  influence  coefficients 
from  an  equation  for  the  induced  velocity.  The  program 
finds  the  strength  of  each  bound  vortex  and  therefore  the 
incremental  circulation  around  the  element  by  solving  the 
set  of  simultaneous  equations  for  wing  shape.  The  program 
can  also  invert  the  problem  and  solve  for  the  wing  shape 
from  a  specified  vortex  distribution. 

The  program  can  then  find  the  forces  and  moments  on  the 
wing,  knowing  the  distribution  of  incremental  circulation. 
This  solution  is  possible  by  using  the  Kut ta - Joukowski 
theorem  which  relates  incremental  circulation  and  effective 
velocity  to  the  force  generated. 


B.  CIRCULATION  MODEL 

This  program  models  the  flow  by  a  series  of  circulation 
elements  distributed  over  the  wing.  There  is  also  a 
circulation  sheet  that  extends  from  the  trailing  edge  to 
inf ini  ty . 

Solution  of  a  set  of  governing  equations  provides  the 
strength  of  each  circulation  element.  The  equations  satisfy 
flow  tangency  on  the  wing  surface  and  two  boundary 
conditions.  Firstly,  circulation  strength  is  zero  along 
the  leading  edge  and  wing  tips.  Secondly,  the  derivative  of 
circulation  with  respect  to  chord-wise  direction  is  zero  at 
the  trailing  edge.  The  derivative  of  circulation  with 
respect  to  chord  is  vorticity.  This  trailing  edge  boundary 
condition  satisfies  the  Kutta  condition  of  zero  vortex 
strength  at  the  trailing  edge.  Zero  vorticity  at  the 
trailing  edge  also  means  that  the  circulation  strength 
behind  the  wing  is  a  function  of  the  span-wise  coordinate 
only  . 

C.  PRESSURE  DISTRIBUTION  MODEL 

The  final  model  of  the  flow  uses  a  series  of  pressure 
difference  elements,  ACp  ,  distributed  over  the  wing.  The 
program  builds  a  set  of  N  equations  in  N  unknowns  from  a 
governing  equation  and  from  the  requirement  for  flow 
tangency.   The  boundary  conditions  are: 

1.   ACp  is  zero  at  the  trailing  edge,  and  along  the  wing 
t ips  . 
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2.  For   a   flat   wing,   the   strength   of   ACp   along   the 
leading  edge  is  unrestrained. 

3.  For   a   wing   of   elliptic   lift   distribution,    the 
strength  of  ACp  along  the  leading  edge  is  zero. 

D.   PROBLEMS  INHERENT  TO  ALL  MODELS 

1 .  Field  and  Control  Point  Placement 

The  placement  of  field  and  control  points  within 
elements  is  important  for  all  models.  Coincident  field  and 
control  points  exhibit  singularities.  Each  model  uses 
different  methods  to  handle  them.  Direct  integration 
eliminates  the  singularity  in  two  models.  The  other 
requires  a  special  arrangement  of  control  and  field  points 
along  the  leading  edge  and  wing  tips. 

2 .  Grid  Size  and  Solution  Speed 

Fine  grids  require  more  time  for  solution  and  larger 
computers.  All  models  compromise  between  data  quantity  and 
speed.  The  circulation  model  is  very  memory  intensive  and 
will  only  run  on  the  mainframe  computer.  The  other  models 
will  run  on  a  desktop  computer,  though  the  time  required 
increases  rapidly  with  finer  grids. 

3 .  Grid  Element  Shapes 

Finite  elements  in  a  non- rec tangular  planform  can  be 
either  trapezoids,  or  rectangles.  Trapezoids  match  the 
planform  exactly;  rectangles  present  a  ragged  leading  and/or 
trailing  edge.  The  form  to  use  is  dependent  on  the  type  of 
model.    The  circulation  model  can  use  a  trapezoid  without 
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difficulty.  Coordinate  transformations  easily  handle  the 
shape  and  at  the  same  time  allow  non-uniform  cosine  spacing 
of  the  elements.  The  Vortex  and  ACp  models  work  best  with 
rectangular  grid  elements.  Therefore  the  planform  is  not  an 
exact  representation  of  the  wing.  This  modification  is 
reasonable  since  the  forces  and  moments  computed  are 
comparable  to  other  models. 

E.    SCOPE  OF  DISCUSSION 

Of  the  three  models  investigated,  only  the  VORTEX 
program  was  completely  successful.  Significant  effort  was 
expended  on  the  CIRCULATION  and  ACp  models  trying  to  make 
them  functional.  Rather  than  expend  a  large  portion  of  this 
thesis  discussing  all  aspects  of  the  two  unsuccessful 
models,  an  overview  of  them  will  be  presented  and  the 
relevant  factors  of  each  will  be  addressed.  Since  the 
VORTEX  program  was  successful,  and  is  completely  functional, 
a  detailed  discussion  of  the  procedure  and  assumptions  is 
provided.  In  keeping  with  this,  the  FORTRAN  source  code  for 
the  VORTEX  model  is  the  only  computer  code  included. 


Ill .  VORTEX  PROGRAM  DEVELOPMENT 

A.   BASIS  OF  THE  VORTEX  MODEL 

The  VORTEX  program  is  an  adaptation  of  a  program  by 
Moran  [Ref  .1]  .  The  program  models  the  flow  by  a  series  of 
horse  shoe  vortices  distributed  over  the  wing.  In  his 
text,  Moran  [Ref.l]  develops  a  simple  program  to  find  the 
strengths  of  a  series  of  horse  shoe  vortices.  His  program 
is  rather  limited  in  that  it  only  works  for  straight  wings 
without  taper  or  sweep.  Additionally,  his  program  only 
finds  the  loading  generated  by  a  flat  wing.  The  final 
limitation  is  his  use  of  uniform  sized  grid  elements,  which 
fails  to  concentrate  grid  elements  near  the  boundaries  of 
the  wing . 

As  with  many  aerodynamic  problems,  the  VORTEX  program 
uses  a  grid  or  mesh  of  finite  sized  elements  distributed 
over  the  surface  of  the  wing.  The  program  divides  the  wing 
planform  (there  is  zero  thickness  modeled  in  this  program) 
into  N  discrete  elements. 

The  program  solves  a  set  of  N  equations  for  N  unknowns. 
The  equations  are  derived  from  an  induced  velocity  equation 
and  the  requirement  for  flow  tangency  on  each  grid  element. 
The  N  unknowns  are  the  strengths  of  the  individual  horse 
shoe  vortices  associated  with  each  element. 


Knowing  the  strengths  of  the  vortex  elements,  the 
program  can  find  the  forces  and  moments  acting  on  the  wing 
as  a  result  of  those  elements. 

B.  MODIFICATIONS  TO  THE  MODEL 

The  VORLAT  program,  developed  by  Moran  [Ref.l],  required 
changes  to  suit  the  goals  of  the  thesis.  Items  requiring 
improvement  were  the  ability  to  handle  sweep  and  taper. 
Other  improvements  were  also  added  to  generate  the  grid 
using  cosine  spacing.  A  module  to  determine  camber  and 
twist,  called  shape,  was  also  incorporated. 

C.  TECHNIQUE  FOR  A  FLAT  WING 

This  is  a  short  description  of  the  VORTEX  program. 
Consult  Moran  [Ref.l]  for  a  description  of  the  "VORLAT 
program,  the  foundation  of  the  VORTEX  program. 

1 .   Down-wash  and  the  Relationship  to  Flow  Tangency 

A  horse  shoe  vortex  is  located  in  the  xy  plane  as 
shown  in  Figure  3-1.  The  vortex  induces  a  velocity  at  any 
point  in  the  xy  plane.  That  velocity  can  be  determined  from 
equation  3.1,  as  follows: 


Ar 

w   (x,y)  =  

V       4n(y-y   ) 


1  + 


\/(x-x  )  +  (y-y  )' 

a  a 

(x-X  ) 
a 


AT 


4n(y-yb) 


1  + 


V(x-x  )  +Qy-jy.) 


(x-x  ) 
a 


3  .  1 
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If  x  =  xa  then, 


Ar  / 
w   (x,y)  =  — 


4n  \  y-ya      y-yt 


3  .  2 


w±j(x,y)  is  the  velocity  induced  at  the  control  point  (x,y). 
AT  is  the  incremental  circulation1,  and  (xa,ya)  and  (xa,yb) 
are  the  coordinates  of  the  corners  of  the  horse  shoe  vortex. 
Moran  derived  this  equation  using  the  Biot-Savart  Law 
[Ref.l].    ^ 


A  Horse  Shoe  Vortex  in  the  Wing  Plane 
Figure  3-1 


Different  authors  use  different  variable  names  for 
Circulation  and  Vorticity.  In  the  Vortex  model,  AT  is  used 
for  incremental  circulation.  In  the  Circulation  model,  T  is 
used  for  Circulation  strength.  Caution  is  urged  when 
relating  one  model  to  another. 
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As  shown  in  Figure  3-2,  the  remote  velocity,  with  a  wing 
at  angle  of  attack  a,  has  components  V^cosCa)     and  VgaSinCa). 


Veo    SirKSO 


Vco  cos«*> 


FLAT    V1NG 


Velocity  Components  on  a  Flat  Wing 
Figure  3-2 

If  the  flow  is  tangent  to  the  wing,   the  down-wash 

o 

velocity^  at  a  point  on  the  wing  must  equal  V^sinCa).  Each 
vortex  on  the  wing  contributes  to  the  down-wash  at  every 
point  on  the  wing.  So,  an  equation  can  be  developed  for  each 
control  point  as  a  function  of  the  strength  of  each  horse 
shoe  vortex.  In  this  thesis,  the  mid  point  of  any  bound 
vortex  is  termed  a  field  point.  Any  point  on  the  wing 
where  flow  tangency  is  evaluated  is  termed  a  control  point. 


^Downwash  is  considered  positive  when  its  direction  is 
downward . 

-^The  contribution  will  be  positive  when  the  induced 
velocity  is  downward,  or  negative  when  the  induced  velocity 
i  s  upward . 
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An  example  of  the  equation  for  the  control  point  (Xp,Yp)     i-s 


V  sina  —  ?    V  w   (x  ,y  )  =  0 

*      —  u    u    p   p 
u 


3  .  3 


Combining  all  the  control  points  results  in  a  linear  set  of 
N  equations  in  N  unknowns. 

2 .   Placement  of  the  Horse  Shoe  Vortex 

Moran  [Ref.l]  and  others  have  discussed  the 
placement  of  a  lumped  vortex  on  a  two-dimensional  airfoil 
with  one  chord-wise  element.  For  the  two-dimensional  case, 
the  vortex  is  placed  at  the  quarter  chord  point  and  flow 
tangency  is  evaluated  at  the  three-quarter  chord  point.  The 
two-dimensional  reasoning  can  be  applied  to  a  wing  of  finite 
span,  assuming  the  wing  is  formed  with  a  single  chord-wise 
e lement . 


o 
o 

DC 
O 


CHORD 


Circulation  Distribution  over  a  Flat  Airfoil 

Figure  3  -  3 
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For  the  flat  airfoil,  the  circulation  distribution 
is  known  to  be  approximately  as  shown  in  Figure  3-3.  The 
centroid  of  this  circulation  distribution  is  at  the  quarter 
chord  point.  As  a  result,  it  is  reasonable  to  lump  the 
total  circulation  at  the  quarter  chord  point.  For  these 
conditions,  the  flow  is  tangent  at  the  three  quarter  chord 
point.  This  arrangement  correctly  computes  the  moment 
coefficient  (Cm)  for  the  two  dimensional  flat  airfoil. 

The  reason  for  other  authors'  placement  of  the 
vortex  at  the  quarter  chord  of  each  grid  element  was 
questioned  during  development  of  the  VORTEX  program.  The 
other  authors  used  lifting  line  theory  with  a  single  chord- 
wise  element  for  their  quarter  chord  vortex  placement.  The 
VORTEX  program  does  not  use  a  single  chord-wise  element. 
The  wing  is  divided  into  a  number  of  individual  elements 
which  are  modeled  with  a  uniform  distribution  of  vorticity 
over  each  element.  The  centroid  is  at  the  center  of  the 
element.  So,  for  the  arrangement  in  Figure  3-4,  the 
incremental  circulation  is  concentrated  at  the  center  of 
each  e lement . 
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Approximation  of  Vorticity  with  Individual  Elements 

Figure  3-4 

With  the  circulation  concentrated  at  the  center  of 
each  element,  the  flow  tangency  point  is  set  at  one  half  an 
element  downstream  of  the  concentrated  circulation  and  falls 
at  the  border  between  elements,  or  at  the  trailing  edge  of 
the  wing. 

The  coefficients  affected  by  the  vortex  placement 
are  C^  and  Xcp .  In  Ref.2,  Xcp/c  for  a  straight  flat  wing, 
aspect  ratio  2,  was  determined  to  be  0.209.  With  the  lumped 
vortex  at  the  mid  chord  point,  the  location  of  XCp/c 
converged  to  a  value  of  approximately  0.234.  The  error 
indicates  the  vortex  is  too  far  aft  on  the  element. 
Relocating  the  lumped  vortex  element  to  the  quarter  chord 
point,  and  evaluating  flow  tangency  at  the  three-quarter 
chord  point,  moved  the  computed  location  of  Xcp/c  forward, 
but  still  not  to  the  true  value  of  0.209. 
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In  Ref.3,  Hough  looked  at  optimum  grid  and  vortex 
arrangement  for  rapid  convergence.  One  of  the  planforms 
was  a  straight  flat  wing  with  an  aspect  ratio  of  2.  Hough 
used  a  conventional  vortex  lattice  arrangement  of  uniform 
size  grid  elements  with  the  lumped  vortex  at  the  quarter 
chord  point.  For  the  aspect  ratio  2  wing,  he  found  that 
the  computed  C^/a  was  larger  than  actual  and  that  the  value 
converged  as  in  Figure  3-5.  It  was  expected  that  the  error 
in  Cl/q:  could  be  reduced  by  decreasing  the  planform  area  by 
some  factor.  He  found  that  convergence  of  Cl/q,  Vortex  Drag 
Factor  (K)  ,  and  Xcp  was  improved  by  insetting  the  tip 
vortex  by  some  fraction  (d)  of  a  grid  element.  This  inset 
improved  convergence  dramatically  when  d  =  1/4.  Though  not 
addressed  by  Moran  [Ref .1] ,  that  is  the  reason  for  insetting 
the  grid  by  d  =  1/4. 

The  VORTEX  program  uses  cosine  spacing  instead  of 
uniform  spacing  so  insetting  the  grid  by  d  =  1/4  was  not  an 
available  option.  Instead,  the  span-wise  grid  layout  is 
developed  with  an  additional  row  of  tip  elements  which  are 
not  used  in  the  computation  of  incremental  circulation  or 
force.  This  reduces  the  planform  area  and  the  improvement 
in   convergence   can   be   seen   in   Figures   3-5   and   3-6. 


Vortex  Drag  Factor,  K  is;   7r«AR»Cni 


CL' 
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Additionally,  the  convergence  of  Xcp  is  improved,  as  seen  in 
Figure  3-7. 
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Figure  3-6 
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Figure  3-7 

3 .  Satisfaction  of  the  Kutta  Condition 

To  satisfy  the  Kutta  condition,  the  vortex  strength 
must  be  zero  at  the  trailing  edge.  While  not  explicitly 
required  by  the  VORTEX  program,  the  resultant  vortex 
strengths,  for  a  flat  wing,  approach  zero  at  the  trailing 
edge.  Therefore,  the  results  seem  to  imply  that  flow 
tangency  at  the  trailing  edge  of  a  flat  wing  is  a  corollary 
of  the  Kutta  condition. 

4 .  Planform  Development 

It  is  necessary  to  develop  wing  geometry  from  aspect 
ratio,  taper  ratio  and  sweep  angle.  Aspect  ratio  is  the 
span  divided  by  the  average  chord.  Taper  ratio  is  the  tip 
chord  divided  by  the  root  chord,  and  sweep  angle  is  the 
angle  between  the  leading  edge  and  the  y  axis.  The  y  axis 
is  perpendicular  to  the  remote  velocity  and  parallel  to  the 
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earth.  The  aircraft  design  courses  taught  at  NPS  use  these 
variables  to  define  the  planform.  Moran  [Ref.l]  uses  a 
fixed  root  chord  equal  to  1.0  and  varies  the  span  to  achieve 
different  aspect  ratios.  The  VORTEX  program  adopts  this 
approach,  which  works  well. 


Remote  Velocity 


> 


b/2 


Planform  and  Variables  used  in  VORTEX  Program 

Figure  3-8 

The  subroutine  SET85  determines  planform  variables 
and  stores  them  in  part  of  the  matrices  WING  and  SECTN  for 
further  use.  The  WING  and  SECTN  matrices  also  contain  final 
results  of  the  program. 

5  .   Grid  Development 

The  program  develops  a  grid  to  model  the  vortex 
system  after  establishing  the  outline  of  the  planform.  The 
grid  elements  are  rectangular.  The  number  of  elements  must 
be  small  to  keep  the  time  required  for  solution  reasonable. 
The  program  concentrates  grid  elements  in  areas  where  the 
rate  of  change  is  most  rapid,  or  where  the  values  are  most 
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important.  Those  areas  are  along  the  boundaries.  An 
excellent  method  for  concentrating  the  points  in  the  areas 
desired  while  minimizing  the  total  number  of  points  is 
through  cosine  spacing.  The  program  uses  a  cosine  function 
to  distribute  span-wise  and  chord-wise  grid  points  in  a  non- 
uniform fashion  after  a  suitable  coordinate  transformation. 
Figure  3-9  contains  a  sample  of  the  layout.  Notice  how  the 
rectangular  elements  model  the  planform. 


Remote  Velocity 


> 


/ 

1       v 

/ 

/ 

\ 

1 

\ 

/ 

/ 

\ 

\ 

\ 

/ 

/ 

^ 

V 

/ 

7^ 

\  1     1 

1  1 

Grid  Model  with  Cosine  Spacing 
Figure  3-9 

6 .   Solving  the  Set  of  Equations 

The  method  used  to  solve  the  set  of  N  linear 
equations  is  arbitrary.  Moran  [Ref .  1]  uses  a  Gaussian 
algorithm  for  solution  of  an  N  by  N  matrix.  The  VORTEX 
program  uses  the  same  algorithm.  Increasing  the  total 
number  of  grid  elements  beyond  100  would  require 
modification  of  the  GAUSS  subroutine. 
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7 .   Finding  the  Forces  on  the  Elements 

The  Kut ta - Joukowski   theorem  states   that   the   force 
per  unit  span  acting  on  an  element  is. 


A/-    >  — * 

—  -  p  V~/x  AT 

Ay  eff 


3.4 


In  this  equation,  Veff  Is  the  local  effective  velocity  at 
the  center  of  the  element.  V eff  is  defined  in  equation  3.5 
and  Figure  3-10.  p  is  the  density  and  AT  is  the  incremental 
circulation  around  the  element.  This  circulation  (AT)  is 
the  same  AT  contained  in  equation  3.1,  and  can  be  termed  the 
incremental  circulation  that  occurs  over  the  element. 
Throughout  this  section,  certain  approximations  and 
dimensional  simplifications  will  be  made.  The  first- is  the 
small  angle  approximation,  where  the  sine  is  approximately 
equal  to  the  angle  in  radians,  and  the  cosine  is 
approximately  equal  to  1.0.  This  approximation  requires 
that  angle  of  attack  for  the  flat  wing  be  small,  generally 
less  than  10  degrees.  The  approximation  also  requires  that 
any  wing  slope  on  the  elliptically  loaded  wing  be  small. 
This  requirement  is  met  by  restricting  the  desired  lift 
coefficient  to  values  less  than  about  0.5.  In  the  process 
of  determining  coefficients  of  lift,  drag  and  moment, 
certain  dimensional  reference  quantities  arise.  These 
quantities  are  density  (p),  remote  velocity  (V^),  planform 
area  (S),  and  average  chord  (c)  .  In  performing  dimensional 
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analysis  an  arbitrarily  value  of  unity  can  be  assigned  to 
thes  e  terms . 

a.   Finding  ACp 

The  distribution  of  ACp  over  the  wing  is 
desired,  so  the  force  on  each  element  must  be  converted  to  a 
dimens ionles s  pressure  difference  coefficient.  ^eff  *-s 
defined  relative  to  the  wing  as; 


V        =  V  cos  a  i    +  V  sin  a  —  w  \  j 


"? 


3  .  5 


Veff 


Components  of  Effective  Velocity  (Veff) 
Figure  3-10 

It  is  also  significant  to  note  that  AT  can  be  written  as; 


A?=Ar? 


3  .  6 
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since  the  bound  portion  of  the  horse  shoe  vortex  is  aligned 
with  the  y  axis.  The  result  of  the  force  cross  product, 
Equa t  ion  3.4,  is ; 


A? 

Ay 


=  P 


-»  -> 

w  -  V   sinal&ri    +Vcosa&Tk 


3  .7 


Setting  p     «=  1  ,  V, 
approximation; 


1,   and  incorporating  the  small  angle 


A£ 

Av 


w  -  a  )AV  i    +   AT  k 


3  .  8 


This  force  is  related  to  ACp  .  ACp  is  a  scalar,  rather  than 
a  vector  quantity,  so  the  program  uses  the  force  component 
normal  to  the  wing  to  compute  ACp  .   That  component  is*; 


AF  =  Ay AT 


3  .  9 


Making   the   force   non-dimensional   gives   a   pressure 
difference  coefficient,  ACp  . 


AC  = 


AF 


P        1 


^pV   AxAy 
2   " 


3  .  10 
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Finally,  after  setting  p    —    1  and  Vffl  =  1,  ACp  for  an  element 


is  ; 


2AF   2Ar 

AC  =  =  

P       AxAy    Ax 


3  .11 


This  equation  provides  the  ACp  for  each  element,  which  can 
be  plotted  versus  chord,  for  each  span-wise  section. 
b.   Finding  the  lift,  drag  and  moment 

Lift,  drag  and  moment  are  also  found  from  the 
vector  force  on  each  element.  Equation  3.7  showed  that  the 
force  on  an  element  has  components  normal  and  tangent  to  the 
wing.  Once  more  setting  p  —  1,  and  V &  =  1,  the  vector 
force  on  an  element  is; 


s>  = 


u>  —  sin  a  lAI'Av  i   +  cos  a  ATAy^ 


3  .  12 


These  forces  are  oriented  relative  to  the  wing  coordinate 
system,  where  i  is  parallel  to  the  wing  and  k  is 
perpendicular  to  the  wing.  Lift  and  drag  are  normal  and 
tangent,  respectively,  to  the  remote  velocity  (Voo).  The 
wing  is  at  an  angle  of  attack  (a)  to  that  remote  velocity. 
Using  this  information,  the  program  transforms  the  forces  to 
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a  coordinate  system  oriented  to  the  remote  velocity.  For  an 
individual  element  the  lift  and  drag  forces  are; 

AL  =  Elemental  Lift  Force  =  (cos  a)  (AyAD  (cos  a)  -  (w  -  sina)(AyAF)  (sin  a) 
AD  =  Elemental  Drag  Force  =  (cosa)(AjyAD  (sina)  +  (w  -  sina)  (AyAn(cosa) 

using  a  small  angle  approximation,  and  discarding 

higher  order  terms,  3.13 

AL  =  Elemental  Lift  Force  -  A>Ar  -  umAjAT  =  AjyAT 
AD  =  Elemental  Drag  Force  —  AyAFa  +  (w  —  a)A_yAr  =  w&y&r 

The  total  lift  and  total  drag  acting  on  the  wing  is  a 
summation  of  the  elemental  lifts  and  drags.  Using  wing 
symme  try ; 


L  =  Total  Lift  =  2        Y        A7. 

semi  -  span  3.14 

D  -  Total  Drag  =  2        V        AD 
semi  —span 


The  moment  about  the  y  axis  for  an  individual 
element,  with  the  accepted  sign  convention,  is; 


AM  =  Elemental  Moment  —  —(normal  wing  force)  (x  position  of  center  of  element) 
AM  =  Elemental  Moment  —  (  — AvAD(x      ) 

cen 
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Total  moment  about  the  y  axis  is  a  summation  of 
the  elemental  moments.   Once  again,  using  wing  symmetry. 


M  =  Total  Moment  -  2  AM  3     16 

semi  —span 


8 .   Finding  the  Coefficients 

Non- dimensionaliz ing  the  total  lift  and  drag  gives, 


C,  = 

-pV2S 
2   °° 


C  = 


3  .  17 
D 


D    »    2 

2   " 

where  S  =  span  x  average  chord  =  (b  x  c)  ,  p  =  1  and  V,,,  -  1. 

Non - dimens ional iz ing   the   total   moment   about"  the  y 
axi s  give s  , 


C,.    = 


M 


Mo       1  s  3.18 

1   9  - 

-pVcS 

where  c  is  the  average  chord. 

The  aerodynamic  center  Xac  is  important.  Xac  is  the 
chord-wise  position  about  which  the  flat  wing  has  zero 
moment,  and  is  found  from. 


—  total  moment  3.19 

00  toto/Zi/i 
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This  equation  assumes  that  total  moment  is  taken  about  the  y 
axis,  or  x  =  0  .  When  the  program  finds  the  moment  on  the 
elliptically  loaded  wing,  it  uses  the  Xac  from  the  flat  wing 
as  the  reference  axis  and  finds  C^ac  . 

D.   TECHNIQUE  FOR  SPECIFIED  LOADING 

1 .   Finding  Camber  and  Twist  from  Loading 

Elliptic  loading  in  a  span-wise  direction  results  in 
minimum  induced  drag,  and  loading  is  proportional  to 
circulation.  At  sufficiently  high  aspect  ratio,  a  flat  wing 
with  elliptic  area  distribution  generates  elliptic  loading 
but  weighs  more  than  a  straight  wing.  Wings  with  straight 
leading  and  trailing  edges  also  cost  less  to  manufacture 
than  elliptic  wings.  So,  the  program  generates  a  span-wise 
elliptic  lift  distribution  by  slightly  twisting  the  wing. 
The  direct  relationship  between  wing  shape  (camber  and 
twist)  and  circulation  distribution  makes  this  generation 
possible  . 
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camber 


twist 


Camber    and    Twist    on    the    Elliptically    Loaded    Wing 

Figure    3-10 

For  a  thin  wing,  elliptic  load  distribution  in  the 
chord-wise  direction  requires  an  approximately  parabolic 
camber  shape.  The  program  specifies  an  elliptic  chord-wise 
load  distribution  and  then  finds  the  wings'  associated 
shape . 

2  .   Specifying  the  Lift  Coefficient 

Specifying  the  form  and  amplitude  of  the  ideal  lift 
distribution  over  the  wing  determines  the  form  and  amplitude 
of  the  wing  shape.  A  load  distribution  of  elliptic  form  is 
imposed  in  this  case  and  the  corresponding  amplitude  is 
fixed  by  the  desired  ideal  wing  lift  coefficient,  C^  • 
3 .    Maximum  Desired  Lift  Coefficient 

The  small  angle  approximation  breaks  down  when  large 
lift  coefficients  are  specified.  To  prevent  this 
occurrence,  the  maximum  desired  lift  coefficient  is 
restricted  to  0.5. 
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4 .  Achieving  the  Desired  Lift  Coefficient 

The  program  scales  a  reference  distribution  of  ACp  , 
which  is  elliptic  in  form,  to  achieve  the  desired  total  lift 
coefficient  for  the  wing. 

5 .  Forces.  Moments  and  Coefficients 

Once  the  ACp  distribution  is  scaled  for  the 
elliptically  loaded  wing,  the  program  finds  forces,  moments 
and  coefficients  in  the  same  manner  as  for  the  flat  wing. 
The  only  difference  is  in  the  coordinate  transformation  used 
to  convert  forces  on  the  wing  to  lift  and  drag.  The 
individual  elements  on  the  elliptically  loaded  wing  are  no 
longer  at  a  uniform  angle  to  the  remote  velocity.  This 
variation  in  angle  must  be  taken  into  account  when 
performing  the  coordinate  transformation  which  gives  lift 
and  drag  forces  on  individual  elements. 

The  program  finds  moment  coefficient  about  the 
aerodynamic  center,"  Cj^ac  ,  from, 
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Mac 


CM      +C, 

Mo  Li 


3  .  20 


where  Xac  is  found  from  the  flat  wing,  and  c  is  the  average 
chord.  Cjvj0  is  the  pitching  moment  of  the  elliptically 
loaded  wing  about  the  y  axis,  and  CLi  is  the  lift 
coefficient  of  the  elliptically  loaded  wing. 
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E.   VALIDATION  OF  RESULTS 

Two  sources  are  used  to  check  the  validity  of  the  VORTEX 
program  results.  Lifting  line  theory  provides  one  source  of 
predicted  lift  and  drag  values  for  straight  high  aspect 
ratio  wings.  A  continuous  loading  method  developed  by  the 
National  Aerospace  Laboratory  of  the  Netherlands  (NLR), 
presented  in  Ref.2,  is  used  as  the  second  source.  The 
method  used  by  NLR  is  a  very  accurate  computational  method 
and  is  considered  to  be  exact  for  this  comparison. 

Kuethe  and  Chow  [Ref .4]  discuss  lifting  line  theory  for 
the  case  of  a  flat,  untapered,  rectangular  wing.  Using 
lifting  line  theory  and  the  VORTEX  program,  lift  and  drag 
coefficients  for  identical  wings  were  computed.  A  range  of 
aspect  rations  from  0.5  through  20  were  selected.  Lifting 
line  theory  is  known  to  be  reliable  at  higher  aspect  ratios, 
but  tends  to  lose  accuracy  at  low  aspect  ratios. 

As  the  graphic  results  in  Figures  3-11  and  3-12 
demonstrate,  CL/a  and  CDi/(a)2  values  obtained  from  lifting 
line  theory  and  the  VORTEX  program  converge  nicely  at  high 
aspect  ratios.  At  lower  aspect  ratio,  lifting  line  and  the 
VORTEX  program  show  a  difference  in  calculated  value.  The 
values  of  CL/a  and  CDi/(a)2  obtained  by  NLR  [Ref.2]  agree 
closely  with  the  VORTEX  program  results  at  aspect  ratio  2. 
The  NLR  results  validate  the  VORTEX  program  results  for 
straight  wings.  Exact  data  from  swept  and  tapered  wings 
were  not  investigated. 
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Comparison    of  VORTEX    Program    results 
with    Lifting    Line   Theory 
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Figure    3-11 

Comparison    of  VORTEX    Program    results 
with    Lifting    Line   Theory 
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IV.  CIRCULATION  MODEL 

A.  BASIS  OF  THE  CIRCULATION  MODEL 

This  program  models  the  flow  with  a  sheet  of  distributed 
circulation  and  a  continuous  trailing  vortex  sheet  behind 
the  wing.  Barna  [Ref .5]  ,  and  also  Milne  -  Thomson  [Ref .6, 
pages  171-177],  discuss  this  conceptually  simple  model.  The 
program  uses  the  wing  shape,  flow  tangency,  boundary 
conditions,  and  the  Biot-Savart  Law  to  develop  a  set  of 
equations.  The  program  then  solves  the  set  of  equations  for 
the  unknown  circulation  strengths  at  all  grid  points  on  the 
wing  . 

B.  TECHNIQUE  FOR  SOLUTION  OF  THE  CIRCULATION  MODEL 
1 .   Induced  Velocity 

The  velocity  induced  at  a  point  by  the  distributed 
circulation  sheet  on  the  wing  and  the  trailing  vortex  sheet 
can  be  treated  as  the  sum  of  the  velocity  induced  by  each. 
The  velocity  induced  at  the  control  point  (x0,y0)  by  field 
point  (x,y),     located  on  the  wing  is. 
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r  is  the  distance  between  the  field  and  control  points,  and 
r  is  the  strength  of  the  circulation  at  the  field  point 
(x,y).  Setting   the   semi-span  length  equal   to  unity,   and 
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varying  the  chord  accommodates  different  aspect  and  taper 
ratios . 


Remote  Velocity 


> 


—  Cn(l-T) 


b/2  =  1 


Cn 


Planform  and  Variables  used  in  Circulation  Model 

Figure  4-1 

Circulation  (T),  and  pressure  difference  coefficient 

(ACD)  are  related  as  follows; 


AC  = 
p 


V 


4.2 


In  equation  4.2,  ACp  is  the  pressure  difference  coefficient 
between  the  upper  and  lower  surfaces,  and  T  is  the 
c  irculat ion . 
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The  velocity  induced  at  control  point  (x0,y0)  by  the 
segment  of  the  trailing  vortex  sheet  attached  to  the  wing  at 
(xt ,y),     is : 


W 


trail 


+  1 


1 


_1     V  4nr  I       (y  —  y  ) 


r—(x„,  —  x  )  i  /  ar„ 

7        oil        T 


[t,dy 


4.  3 


r  is  the  distance  between  the  points,  and  Tt  is  the  strength 

of  the  circulation  at  the  field  point  (X(-,y)     on  the  trailing 
e  dge  . 

When   the   remote   velocity   is   unity,  the   induced 

velocity  is  the  same  as  the  slope  of  the  wing  at  the  point 
in  question,  and  the  total  slope  is: 
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The   program   uses   these   equations   to   satisfy   flow 
tangency  on  the  wing  with  a  known  shape  and  solves  for  the 
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unknown   circulation   strengths.     In   matrix   notation   the 
equations  have  the  following  form. 


ar 


dx 


+ 


B 


ay 


dz 
dx 


4.  7 


2  .   Boundary  Conditions 

The  requirement  for  flow  tangency  at  each  element  on 
the  wing  is  not  enough  to  find  the  circulation  which 
satisfies  the  Kutta  condition;  the  program  must  also  enforce 
certain  boundary  conditions. 

The  flat  wing,  which  produces  the  additional  lift, 
has  the  following  boundary  conditions: 

1.  T  =  0  along  the  leading  edge  and  the  wing  tips. 

2.  dT/dx    -  0  along  the  trailing  edge. 

The  elliptically  loaded  wing,  which  produces  the 
ideal  lift,  has  the  following  boundary  conditions: 

1.  r  -  0  along  the  leading  edge,  and  wing  tips. 

2.  dT/dx    =  0  along  both  the  leading  and  trailing  edges. 

Using  a  finite  difference  scheme,  it  is  possible  to 
approximate  the  partial  derivatives  of  T  (dT/dx,  &  dT/dy) 
from  the  values  of  T  at  neighboring  elements  and  satisfy 
these  boundary  conditions.  The  partial  derivative  of  T  with 
respect  to  x,  namely  (dT/dx),  approaches  infinity  at  the 
leading  edge  of  the  flat  wing.  This  singularity  makes 
finite  differencing  near  the  leading  edge  difficult. 
Cosine  spacing  is  used  to  minimize  that  problem  and  also  to 
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spacing  is  used  to  minimize  that  problem  and  also  to 
concentrate  elements  near  the  edges  of  the  wing.  The  x 
coordinate  in  the  chord-wise  direction  can  be  expressed  in 
an  alternate  form  by  <f>  ,  where: 


(}>  =  cos 
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4.  8 


The  additional  variables  are  defined  in  equation  4.24  and 
Figure  4-3.  The  y  coordinate  in  the  span-wise  direction 
becomes  8,     where: 


8  =  cos     (y) 


4.  9 


After  this  transformation 
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4.10 


From  the  transformation,  it  is  possible  to  show  that  dT/d<f> 
has  a  finite  value  at  the  leading  edge  even  when  dT  /  dx  is 
infinite  . 


36 


0 

Circulation  vs  Phi 
Figure  4-2 

Using  Figure  4-2  as  a  guide,   the  program  estimates 

the   value   of  dT/d<j>      from   the   values   of  T      on  neighboring 

elements.    The   following  shows   the  procedure   for   finding 

8T/d4>     at   the   leading  edge   element  as  a  function  of  T     on 

neighboring  elements.   Other  points  are  derived  in  a  similar 

fashion . 
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For  the  first  element 


For  the  second  element, 
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For    the    third    element, 
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dr  \  A$  /  A<J>  V 
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when  4.16  is  combined  with  4.15,  it  becomes; 
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The  uncertainty  in  this  approximation  is  of  order  A<£4  .   For 
interior  elements, 
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and  for  the  trailing  edge  element, 
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Similar   relations   hold   in   the   span-wise   direction.     In 
matrix  form: 
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3 .   Sweep  and  Taper 

The  program  must  model  swept  and  tapered  wings.  It 
makes  a  series  of  coordinate  transformations  to  accomplish 
this.  The  final  result  is  an  equation  for  the  slope  due  to 
the  circulation  on  the  wing. 
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and  for  the  contribution  of  the  trailing  vortex  filament, 
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where 


x  =  (x  -  fry) 


X    =  (x    -  Pj  ) 
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\  =  Lan(A) ,  tangent  of  leading  edge  sweep 
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4 .   Field  and  Control  Point  Placement 

Field  and  control  points  are  located  at  the  center 
of  the  grid  elements.  When  each  grid  element  contains  both 
a  field  and  control  point,  the  circulation  that  results  from 
solving   the   equations   oscillates   wildly   and   bears   no 
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resemblance  to  the  expected  solution.  After  considerable 
probing,  a  modified  arrangement  was  discovered  which 
generates  a  reasonable  function.  Distribution  of  a  "curb" 
of  elements  containing  only  field  points  along  the  tips  and 
leading  edge  gives  a  satisfactory  solution.  The  remainder 
of  the  grid  elements  contain  field  and  control  points.  The 
final  configuration  is. 


Renote   Velocity 
> 


Field   and   Control   Points 


Field   Points   Only 


Modified  Grid  Geometry 
Figure  4-3 

This  arrangement  gives  up  some  degrees  of  freedom,  but  that 

can  be  compensated  by  adding  an  extra  row  of  grid  elements 

along  the  chord  and  an  extra  column  of  grid  elements  along 

the  span. 

5 .   Influence  Coefficient  Matrices 

Equations  4.22  and  4.23  can  be  divided  into  multiple 

integrals.     The   partial   derivatives  (8T/d4>      &  dT/86)      are 

evaluated  at  the  centers  of  the  elements,  and  the  integrals 

in  4.22  and  4.23  can  be  represented  as  a  summation  over  the 

elements.     Using   these   assumptions   and   equations   4.19 
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through  4.24  it  is  possible  to  write  the  set  of  equations  in 
matrix  form . 

The  program  builds  the  final  matrix  of  influence 
coefficients  from  a  number  of  subsidiary  matrices.  When 
equation  4.20  and  4.21  are  combined  with  equation  4.22  and 
the  cosine  coordinate  transformation  we  get. 


M 
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4.25 


Once   [M]   is   generated,   it   is   possible   to   obtain   the 
circulation  solution  vector  {T). 

6 .  Solving  the  Set  of  Equations 

The  program  uses  the  LEQIF  subroutine  from  the  IMSL 
FORTRAN  library  as  a  linear  equation  solver. 

7 .  Finding  the  Lift.  Drag,  and  Moment 

From  the  distribution  of  circulation  strength  it  is 
possible  to  find  the  lift  and  drag.  Barna  [Ref.5]  and 
Mi lne -  Thorns  on  [Ref.6]  derive  the  equations  for  lift  and 
drag,  as  a  function  of  circulation  and  down-wash.  The 
program  uses  these  equations  to  get  lift  and  drag  from  the 
circulation  on  the  trailing  edge  elements.  The  lift  on  a 
span-wise  section  which  lies  between  y  and  (y+Ay)  is, 


AL  =  pVFAv 


4.  26 
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and  the  induced  drag  on  a  section  is. 

AD  = pVVwAy  ^  .  27 

T  is  the  circulation  along  the  trailing  edge  and  w  is  the 
down-wash  at  the  trailing  edge.  The  procedure  for  finding 
moment  and  center  of  pressure  is  not  as  simple.  The  program 
would  need  to  resolve  forces  on  the  individual  elements  to 
find  moment  coefficient  and  center  of  pressure.  The 
CIRCULATION  program  was  abandoned  before  Fortran  code  to 
generate  moment  coefficient  and  center  of  pressure  could  be 
wr  i  t ten . 

C.   PROBLEMS  ENCOUNTERED 

1 .   Coincident  vs  Non- co inc ident  Points 

The  manner  used  to  handle  singularities  is  one  of 
the  biggest  problems  faced  when  using  finite  elements  to 
model  continuous  functions.  The  circulation  model  behaves 
very  well  when  the  field  and  control  elements  are  not 
co  inc  ident  . 

Field  elements  induce  velocities  which  increase  very 
rapidly  as  the  distance  to  the  control  point  decreases. 
When  the  field  and  control  points  are  in  the  same  element  it 
appears  the  induced  velocity  is  indeterminate.  However, 
Milne -Thomson  [Ref.6]  shows  that  an  element  induces  zero 
velocity  on  itself.  This  result  of  the  singularity  is 
important  and  requires  special  attention. 
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2 .   Field  and  Control  Points  in  the  Wing  Interior 

The  direction  of  velocity  induced  at  a  control  point 
by  its  neighboring  field  points  is  significant.  Using  a 
simple  3X3  element  grid  near  the  center  of  the  wing  as  an 
example,  Figure  4-4  shows  the  direction  of  the  induced 
velocity  from  8  field  points.  All  eight  field  points  have 
positive  circulation  and  surround  the  center  control  point. 
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Interior  Control  Point 
Figure  4-4 

The  three  elements  upstream  from  the  central  control  point 

induce  a  down-wash  on  the  control  point.    The  other  five 

elements  induce  an  up-wash  on  the  control  point.    The  net 

velocity  is  a  sum  of  all  eight  elements.    Recall  that  the 

velocity  induced  is  dependant  on  dT/dx,     not  on  T  itself1. 

This  net  velocity  will  be  down-wash  if  dT/dx     in  the  three 

upstream  elements  is  sufficiently  greater  than  in  the  other 

five  elements . 


The  derivative  of  circulation  is  proportional  to  ACp. 
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3 .   Field  and  Control  Points  at  the  Leading  Edge 

A  problem  becomes  apparent  when  checking  the  induced 
velocity  at  an  element  on  the  leading  edge.  For  a  2X3  grid 
on  the  leading  edge,  Figure  4-5  shows  the  direction  of 
velocity  induced  by  each  of  the  neighboring  field  points  on 
the  central  control  point. 
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control   point 


Leading  Edge  Control  Point 
Figure  4-5 

For   this   situation,   all   the   induced  velocity   is   up-wash, 

assuming  positive  circulation.   However,  the  wing  requires  a 

net  down-wash  for  flow  tangency.   Using  this  arrangement  of 

points,   the   solution   vector   oscillates   wildly   between 

positive  and  negative  values.   This  oscillation  is  a  result 

of  the  system  trying  to  generate  down-wash  at  the  leading 

edge  e 1 ement  s  . 

The  program  uses  a  simple  solution  to  this  problem. 

By  ignoring  the  requirement  for  flow  tangency  on  the  row  of 

elements  along  the  tips  and  leading  edge,  the  program  gives 

a  net  down-wash  at  all  other  elements.    This  simplification 

results   in  a  well  behaved,   positive   circulation  over   the 

ent  ire  wing . 
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It  appears  the  set  of  equations  has  more  unknowns 
than  equations,  but  boundary  conditions  provide  the 
additional  constraints. 

4 .   Model  Complexity 

Developing  a  computer  program  usable  as  a  teaching 
tool  for  basic  aerodynamics  is  the  principal  goal  of  the 
thesis.  The  CIRCULATION  program  does  not  fully  achieve  that 
goal.  The  concept  is  relatively  understandable,  but  is 
difficult  to  execute.  The  complex  method  needed  to  build 
the  influence  coefficient  matrix  reduces  its  usefulness  as  a 
teaching  tool.  The  student  needs  a  simpler  tool.  The 
VORTEX  program  discussed  in  Section  III  is  that  tool. 
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V.  PRESSURE  DIFFERENCE  MODEL 

A.    BASIS  OF  THE  ACp  MODEL 

Like  the  vortex  and  circulation  models,  the  pressure 
difference  model  (ACp)  also  uses  a  governing  integral 
equation  as  its  foundation.  The  integral  equation  relates 
wing  slope  at  a  specific  point  on  the  wing  to  the 
distribution  of  pressure  difference  over  the  whole  wing. 
Dividing  the  wing  into  N  grid  elements,  with  known  slopes, 
the  program  solves  for  N  pressure  differences  (ACp). 

The  requirements  of  the  Kutta  condition  are  not 
explicitly  satisfied  by  the  equation.  Instead,  additional 
constraints  are  necessary  at  the  wing  tips  and  trailing 
edge.  The  simplest  form  of  satisfaction  is  to  require  that 
ACp  be  zero  at  those  grid  elements.  That  constraint  is  more 
severe  than  necessary.  ACp  must  be  zero  at  the  edge  of  the 
element  but  the  value  at  the  center  of  the  edge  elements  can 
be  finite  and  is  approximated  by  fitting  a  polynomial 
through  the  edge  and  neighboring  control  points.  This 
approach  is  similar  to  that  used  in  the  CIRCULATION  model. 
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B.   TECHNIQUE  FOR  SOLUTION  OF  THE  ACp  MODEL 
1 .   The  Singularity 

The    integral   equation   governing 
slope/pressure  difference  relationship  is, 


the        wing 
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where    the    slope, 
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When  the  field  and  control  points  are  coincident  (r=0)  or 
share  the  same  y  value  (y=~yl),  there  is  a  strong 
singularity.  It  is  possible  to  evaluate  this  integral  and 
eliminate  the  singularity. 
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2 .   Evaluation  of  the  Integral 

Using  the  mirror  image  element  on  the  opposite  semi- 
span,  it  is  possible  to  reduce  the  integral  to, 


Z°(x,y) 


S/2 


k0°(x-xl,y-yl)&Cp(xl,yl)dxldy1 


5.4 


where , 


and 


k    =  k.+  k 


5  .  5 


k   = 


1    1 


8n  Cv+jy,)2 


l  + 


(x-ij) 


y/(x-xlf  +  (y+yJ 


5.  6 


The  singularity  was  eliminated  by  evaluating  the  integral 
Using  Figure  5-1  as  a  guide,  the  result  is. 
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XI 

Variables  used  for  a  Grid  Element 
Figure  5-1 
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k=7Z^%+^o  +  ^u  +  ^2) 


8n 


Q^ix-xJ  +  S 


Q2  =  (x-x,)-6 


Q3  =  (y  -  y,)  +  c 


Q4  =  (y-yj)-e 


Q5=VQ2  +  Q2 


*2    ,    ^.2 


q6  =  V<^  +  q; 


^2    ,    ^,2 


Q7  =  V^  +  <j; 


q8=vq;+«; 


«10  = 


4cS 


Q, 


Q, 


«n  = 


-Q4 


Q, 


In 


Q* 


+  Qt 


Q. 


<?12  = 


+  Q. 


<?3<?4 


Q, 


+ 


Qr 


«3«4 


-  1      /n 


■Zn 


<?3       +^5 


Q, 


+  Q, 


and  k" ,     the  mirror  image  is; 


*°=^(Q9+Qio  +  Sn  +  Q12) 
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where , 


Q3  =  (y  +  >-,)  +  c 

Q4  =  (y  +  yj)-c 

and  all  other  values  are  the  same  as  for  k.        The  end  result 
is  , 


k  =  k  +  k 


and 


Z\x,y)=   V  rAC^Xjj,) 
S/2 


In  matrix  form,  the  equations  can  be  written. 


ACp 


A   computer   program   will   rapidly   calculate   all   of   these 
value  s . 

3 .   The  Kutta  Condition 

Solution  of  the  equations  does  not  guarantee  a 
result  which  satisfies  the  Kutta  condition  .  To  enforce  the 
Kutta  condition,  the  program  sets  ACp  at  elements  along  the 


1  It  should  be  noted  that  the  requirement  for  flow 
tangency  at  the  elements  near  the  trailing  edge  was 
sufficient  to  satisfy  the  Kutta  condition  in  the  VORTEX 
program.  However  that  result  was  discovered  after  the 
PRESSURE  program  was  abandoned  and  no  effort  was  made  to 
extend  that  reasoning  to  the  PRESSURE  program. 
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wing  tip  and  trailing  edge   (the  nth  elements)  equal  to  a 

fraction  of  the  value  of  the  [n-l]th  elements.  The  program 

enforces  the  condition  by  fitting  a  polynomial  through  the 
two  points.  See  Figure  5-2. 


REMDTE  VELOCITY 


^> 


N 

N-l 

N  and  N-lth  Elements 
Figure  5-2 

4 .   Grid  Spacing 

The  program  uses  a  non-uniform  grid  spacing,  based 

on  a  cosine  function.    The  elements  near  the  wing  tip  are 

small.    This  spacing  has  the  advantages  noted  earlier. 

C.   PROBLEMS  ENCOUNTERED 

In  matrix  form,  the  model  has  the  form. 


ACp} 


5  .7 


k° '  is  the  influence  coefficient  matrix,  and  z"  is  the  wing 
slope.    The  model  specifies  some  values  of  ACp  on  the  left 
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side  of  the  equality  and  some  values  of  wing  slope  (z°)  on 
the  right  side  of  the  equality.  This  complexity  prevents 
use  of  the  normal  matrix  solvers.  As  a  result,  the  solution 
requires  a  complex  matrix  manipulation  which  is  not  readily 
understandable,  or  desirable  for  a  teaching  tool. 

D.   CONCLUSIONS 

The  evaluation  of  the  integral  over  the  element  has 
definite  advantages  for  coincident  field  and  control  points. 
The  evaluation  eliminates  a  strong  singularity  and  should 
give  a  result  using  finite  elements  that  is  very  close  to 
the  actual  property. 

The  Kutta  condition  and  method  used  to  enforce 
satisfaction  is  not  optimal.   Further  studies  are  needed. 
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B.    INTRODUCTION 

LCDR  Chris  L.  HOLM  wrote  this  program,  as  a  tool  for  use 
in  AE2035.  The  program  is  partial  satisfaction  of  the 
requirements  for  the  degree  Aeronautical  Engineer  from  the 
Naval  Postgraduate  School.  The  objective  of  the  program  is 
to  provide  a  simple  computer  simulation  of  flow  over  a  thin 
wing  with  low  aspect  ratio.  High  aspect  ratio  wings  are 
better  served  by  a  lifting  line  model.  Many  advanced 
programs  using  elaborate  methods  to  model  viscous  effects, 
compressibility,  boundary  layer  growth,  and  shocks,  are 
available.  None  of  them  gives  a  good  introduction  to  the 
field  of  computational  aerodynamics.  This  program  should 
fill  that  need  at  the  Naval  Postgraduate  School. 

If  you  wish  to  skip  the  theory  behind  the  program,  go  to 
section  I  for  instruction  on  running  the  program. 

There  are  a  number  of  ways  to  model  the  source  of 
forces  acting  on  a  body  surrounded  by  fluid  in  motion. 
These  include  potential  functions,  vortex  distributions, 
circulation  distributions,  and  pressure  differential 
distributions.  Each  model  is  related  to  the  other,  and  each 
has  advantages  and  disadvantages.  This  program  uses  a  set 
of  horse  shoe  vortices  to  approximate  the  flow  over  a  wing. 
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C.    SYSTEM  REQUIREMENTS 

System  type 

Memory 

Math  Co-processor  (8087) 

Graphics  card 


Graphics  program 


IBM    PC    or    compatible 
256K,     minimum 
Recommended,  not 

required . 

Recommended,  not 

required . 

X,        Y       graphing       program 
like  GRAPHER,  is 

recommended,  not 

required . 

D.   CONCEPT  DESCRIPTION 

1 .   The  Horse  Shoe  Vortex 

This  program,  which  models  the  flow  by  a  series  of 
horse  shoe  vortices,  distributed  over  the  wing,  is  an 
adaptation  of  the  VORLAT  program  by  Moran  [Ref.l].  He 
develops  a  program  to  find  the  strengths  of  a  series  of 
horse  shoe  vortices  associated  with  a  flat  rectangular  wing. 
The  VORLAT  program  has  its  foundation  in  two-dimensional 
airfoil  theory,  and  Moran  [Ref.l]  adapted  the  theory  to 
wings  of  finite  aspect  ratio. 

The  Users  Manual  will  present  a  short  description  of 
the  VORTEX  program.  For  coverage  of  the  VORLAT  program,  the 
foundation  of  the  VORTEX  program,  consult  Moran  [Ref.l]. 
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a.   Down-wash  and  the  Relationship  to  Flow  Tangency 

The  velocity  induced  at  a  point  in  the  xy  plane 
by  a  horse  shoe  vortex,  also  located  in  the  xy  plane,  can  be 
determined  from: 


Ar 

w  ix,y)  =   

U       An(y-y 


1  + 


V(i-i  ?+(y-y  ) 


(x-x  ) 
a 


AT 


4n(y-yb) 


1  + 


VU-xf+ly-yf 


(x-x  ) 
a 


7  .  1 


if  (x    -  xa)  , 


AT  / 

w   Uj)  =  — 


4n  \y-ya      y-yb 


7  .  2 


w    is  the  velocity  induced  at  (x,y),     called  a  control  point. 
AT  is  the  strength  of  the  horse  shoe  vortex,  and  (xa,ya)     and 
(xa>yb)      are   t^ie   corners   of   the   horse   shoe   vortex.     The 
center  of  the  horse  shoe  vortex  is  a  field  point. 

The  velocity  at  a  point  on  a  flat  wing,  which 
is  at  an  angle  of  attack  a,  has  components  V  ^cos  (a)  and 
VggS  in  (a)     as  indicated  in  Figure  7-1. 


V„  sinCd? 


V,^  cos(cO 


FLAT    V1NG 


Velocity  Components  on  a  Flat  Wing 
Figure  7-1 
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In  order  to  ensure  the  flow  is  tangent  to  the 
wing,  the  induced  velocity  or  down-wash1  at  that  point  on 
the  wing  must  also  be  VajSinCa).  Each  vortex  on  the  wing 

contributes  to  the  down-wash^  at  every  control  point  on  the 
wing.  So,  an  equation  for  each  control  point,  as  a  function 
of  all  the  field  points  can  be  written.  An  example  of  the 
equation  for  the  control  point  (xp,yp)     is. 


V  sina   -  ^  Ar  w   (x  ,y  )  =  0  -,  Q 

°°      —   y  V    pjp  7  .  3 

y 


When  the  control  points  are  combined,  the  result  is  a  set  of 
N  equations  in  N  unknowns. 

b.   Placement  of  the  Horse  Shoe  Vortex 

The  bound  portion  of  the  vortex  (that  portion 
perpendicular  to  the  onset  flow)  is  placed  at  the  1/4  chord 
point  of  each  grid  element  on  the  wing.  Flow  tangency  is 
evaluated  at  the  3/4  chord  point  of  each  grid  element. 


Down-wash  is  considered  positive  when  its  direction  is 
downward . 

^The  contribution  will  be  positive  when  the  induced 
velocity  is  downward,  or  negative  when  the  induced  velocity 
i  s  upward . 
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c.   Satisfaction  of  the  Kutta  Condition 

To  satisfy  the  Kutta  condition,  the  vortex 
strength  must  go  to  zero  at  the  trailing  edge.  For  wings 
and  airfoils  of  zero  thickness,  flow  tangency  enforced  at 
the  trailing  edge  appears  to  be  a  corollary  of  the  Kutta 
condi  t  ion . 

2 .   Developing  the  Wing 

a.   Planform  Geometry 

It  is  necessary  to  develop  wing  geometry  from 
aspect  ratio,  taper  ratio  and  sweep  angle.  Aspect  ratio  is 
the  span  divided  by  the  average  chord.  Taper  ratio  is  the 
tip  chord  divided  by  the  root  chord,  and  sweep  angle  is  the 
angle  between  the  leading  edge  and  the  y  axis.  Moran 
[Ref.l]  uses  a  fixed  root  chord  and  varies  the  span  as 
necessary  to  achieve  the  required  aspect  ratio  and  taper 
ratio.  The  VORTEX  program  uses  this  approach,  which  works 
well. 
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Renoie    Velocity 


> 


b/2 


Planform  and  Variables  used  in  VORTEX  Program 

Figure  7-2 

Aspect   ratio,   AR ,   is  (b^/S).  Span  (b)       is   the 

distance  from  wing  tip  to  wing  tip,  and  area  (S)  is  the 
total  planform  area  of  the  wing. 

Taper  ratio,  A,  is  the  ratio  of  tip  chord  to  root 
chord . 


A  = 


tip  chord         ct 
root  chord        c 


7  .4 


Sweep  angle  delta,  A,  is  the  angle  that  the  leading 
edge  makes  with  a  line  drawn  perpendicular  to  the  root 
chord  . 

b.   Grid  Geometry 

Using  aspect  ratio,  taper  ratio,  and  sweep 
angle,  the  subroutine  SET85  develops  a  grid.  The  program 
stores  the  coordinates  of  all  necessary  points  in  the 
matrices  WING  and  SECTN.  Results  of  the  program  are  also 
stored  in  WING  and  SECTN  which  are  written  to  data  files 
when  the  program  finishes. 
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In  a  manner  analogous  to  that  used  by  Hough  [Ref.3], 
the  tip  vortices  are  inset  to  improve  accuracy  of  the 
results.   The  inset  distance  is  one  element  wide. 

The  model  keeps  matrix  size  and  the  number  of  grid 
elements  small  to  keep  the  time  necessary  for  solution 
within  reason.  The  model  concentrates  grid  elements  in 
areas  where  the  rate  of  change  is  rapid,  or  where  the  values 
are  most  important.  Those  areas  are  along  the  wing 
boundaries.  An  excellent  method  for  concentrating  the 
points,  while  keeping  the  total  number  of  points  at  a 
minimum,  is  through  cosine  spacing.  Figure  7-3  contains  a 
sample  of  the  layout.  Notice  that  the  rectangular  elements 
do  not  model  the  planform  exactly,  but  in  the  limit  the 
difference  will  be  negligible. 
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Grid  Model  with  Cosine  Spacing 
Figure  7-3 
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3 .   Matrices  Used  in  VORTEX  Program 

The  program  uses  two  identical  matrices  of 
influence  coefficients,  [A]  and  [ AA ] .  [A]  is  used  to  solve 
for  the  AT  of  the  flat  wing,  and  [AA]  is  used  to  solve  for 
the  wing  shape  of  the  elliptically  loaded  wing. 


-  i 


dx 


—     =     AT         (Flat  Wing) 


7  .5 


A  A 


AT    = 


dz  | 

—         (Elliptically  Loaded  Wing) 
dx  J 


7.  6 


The  SECTN  matrix  contains  values  for  the  wing 
sections,  (chord,  and  width)  as  well  as  some  of  the  final 
results.  Each  row  of  the  SECTN  matrix  corresponds  to  a 
wing  section  in  the  semi-span.  The  Output  Variables 
section  (H.l.a)  describes  the  columns. 

The  WING  matrix  contains  wing  coordinates  required 
by  the  program  as  well  as  some  of  the  final  results.  Each 
row  of  the  matrix  corresponds  to  an  element  in  the  wing 
semi-span.  The  Output  Variables  section  (H.l.b)  describes 
the  c  o lumns . 


■^The   method   used   to   develop   wing   shape   will   be 
discussed  in  section  VII. E. 


66 


4 .  Solution  of  the  Set  of  Equations 

The  subroutine  GAUSS  solves  the  N  equations  in  N 
unknowns.  The  subroutine  destroys  [A]  in  the  process  but 
returns  the  solution  vector  (AT)  in  place  of  the  right  hand 
s  ides . 

5 .  Finding  the  Forces  on  the  Elements 

The  Kutta - Joukowski  theorem  states  that  the  force 
per  unit  span  acting  on  an  element  is. 


AF    >  — > 

Ay  eff 


7  .7 


In  this  equation,  V eff  i-s  *-he  local  effective  velocity  at 
the  center  of  the  element.  ^eff  ^s  defined  in  equation  7.8 
and  Figure  7-4.  p  is  the  density  and  AT  is  the  incremental 
circulation  around  the  element.  This  circulation  (AT)  is 
the  same  AT  contained  in  equation  7.1.  AT  can  be  termed  the 
incremental  circulation  that  occurs  over  the  element. 
Throughout  this  section,  certain  approximations  and 
dimensional  simplifications  will  be  made.  The  first  is  the 
small  angle  approximation,  where  the  sine  is  approximately 
equal  to  the  angle  in  radians,  and  the  cosine  is 
approximately  equal  to  1.0.  This  approximation  requires 
that  angle  of  attack  for  the  flat  wing  be  small,  generally 
less  than  10  degrees.  The  approximation  also  requires  that 
any  wing  slope  on  the  elliptically  loaded  wing  be  small. 
This   requirement   is   met   by   restricting   the   desired   lift 
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coefficient  to  values  less  than  about  0.5.  In  the  process 
of  determining  coefficients  of  lift,  drag  and  moment, 
certain  dimensional  reference  quantities  arise.  These 
quantities  are  density  (p),  remote  velocity  (Veo),  planform 
area  (S),  and  average  chord  (c).  In  performing  dimensional 
analysis  an  arbitrarily  unit  value  can  be  assigned  to  these 
terms . 

a.   Finding  ACp 

The  distribution  of  ACp  over  the  wing  is 
desired,  so  the  force  on  each  element  must  be  converted  to  a 
dimens ionles s  pressure  difference  coefficient.  V  eff  ^s 
de  fined  as . 


>  *        /  \  * 

V  ,,    =  V  cos  a  i    +      V  sin  a  —  w     j 

eff  CD  \       a.  / 


1  .  8 


Veff 


Components  of  Effective  Velocity  (Veff) 
Figure  7-4 
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It  is  also  significant  to  note  that  AT  can  be  written  as, 


Ar  -  AFj 


7.9 


since  the  bound  portion  of  the  horse  shoe  vortex  is  aligned 
with  the  y  axis.  The  result  of  the  force  cross  product, 
Equat  ion  7.7,  is . 


A? 


w  —  V   sin  a  )AF  i    +  V  cos  a  A  T  k 


7.  10 


Setting  p     =  1,  Va,     =  1,   and  incorporating  the  small  angle 
approximation . 


A/-' 

a7 


w  -  a   AH  +ATi 


7.11 


This  force  is  related  to  ACp  .  ACp  is  a  scalar,  rather  than 
a  vector  quantity,  so  the  program  uses  the  force  component 
normal  to  the  wing  to  compute  ACp .   That  component  is. 


AF  =  AyAr 


7.12 


Making        the        force        non-dimensional        gives        a        pressure 
difference    coefficient,     ACp . 


AC   = 
p 


AF 


p^AxAy 


7.13 
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Finally,  after  setting  p    =  1  and  V^  =  1,  ACp  for  an  element 
is  . 


AC 


2  AF   2 Ar 


7.14 


AxA,y    Ax 


This  provides  the  ACp  for  each  element,  which  can  be  plotted 
versus  chord,  for  each  span-wise  section. 

b.   Finding  the  Lift,  Drag  and  Moment 

Lift,  drag  and  moment  are  also  found  from  the 
vector  force  on  each  element.  Equation  7.10  showed  that  the 
force  on  an  element  has  components  normal  and  tangent  to  the 
wing.  Once  more  setting  p  =  1,  and  ]/„  —  1,  the  vector 
force  on  an  element  is. 


AF  = 


w  —  smalAFA^j    +  cosaATAv  k 


1  .15 


These  forces  are  oriented  relative  to  the  wing  coordinate 
system,  where  i  is  parallel  to  the  wing  and  k  is 
perpendicular  to  the  wing.  Lift  and  drag  are  normal  and 
tangent,   respectively,   to   the   remote  velocity  (V«>)-  The 

wing  is  at  an  angle  of  attack  («)  to  that  remote  velocity. 
Using  this  information,  the  program  transforms  the  forces  to 
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wing  tip  and  trailing  edge   (the  nth  elements)  equal  to  a 

fraction  of  the  value  of  the  [n-l]th  elements.  The  program 

enforces  the  condition  by  fitting  a  polynomial  through  the 
two  points.  See  Figure  5-2. 


REMOTE  VELDCITY 


-> 


N 

N-l 

N  and  N-lth  Elements 
Figure  5-2 

4 .   Grid  Spacing 

The  program  uses  a  non-uniform  grid  spacing,  based 

on  a  cosine  function.    The  elements  near  the  wing  tip  are 

small.    This  spacing  has  the  advantages  noted  earlier. 


C.   PROBLEMS  ENCOUNTERED 

In  matrix  form,  the  model  has  the  form 

AC  J  =  \z 


'' 


5.  7 


k° °     is  the  influence  coefficient  matrix,  and  z"     is  the  wing 
slope.    The  model  specifies  some  values  of  ACp  on  the  left 
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side  of  the  equality  and  some  values  of  wing  slope  (z°)  on 
the  right  side  of  the  equality.  This  complexity  prevents 
use  of  the  normal  matrix  solvers.  As  a  result,  the  solution 
requires  a  complex  matrix  manipulation  which  is  not  readily 
understandable,  or  desirable  for  a  teaching  tool. 

D.   CONCLUSIONS 

The  evaluation  of  the  integral  over  the  element  has 
definite  advantages  for  coincident  field  and  control  points. 
The  evaluation  eliminates  a  strong  singularity  and  should 
give  a  result  using  finite  elements  that  is  very  close  to 
the  actual  property. 

The  Kutta  condition  and  method  used  to  enforce 
satisfaction  is  not  optimal.   Further  studies  are  needed. 
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VI.  VORTEX  PROGRAM  FLOW  CHART 
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B.    INTRODUCTION 

LCDR  Chris  L.  HOLM  wrote  this  program,  as  a  tool  for  use 
in  AE2035.  The  program  is  partial  satisfaction  of  the 
requirements  for  the  degree  Aeronautical  Engineer  from  the 
Naval  Postgraduate  School.  The  objective  of  the  program  is 
to  provide  a  simple  computer  simulation  of  flow  over  a  thin 
wing  with  low  aspect  ratio.  High  aspect  ratio  wings  are 
better  served  by  a  lifting  line  model.  Many  advanced 
programs  using  elaborate  methods  to  model  viscous  effects, 
compressibility,  boundary  layer  growth,  and  shocks,  are 
available.  None  of  them  gives  a  good  introduction  to  the 
field  of  computational  aerodynamics.  This  program  should 
fill  that  need  at  the  Naval  Postgraduate  School. 

If  you  wish  to  skip  the  theory  behind  the  program,  go  to 
section  I  for  instruction  on  running  the  program. 

There  are  a  number  of  ways  to  model  the  source  of 
forces  acting  on  a  body  surrounded  by  fluid  in  motion. 
These  include  potential  functions,  vortex  distributions, 
circulation  distributions,  and  pressure  differential 
distributions.  Each  model  is  related  to  the  other,  and  each 
has  advantages  and  disadvantages.  This  program  uses  a  set 
of  horse  shoe  vortices  to  approximate  the  flow  over  a  wing. 
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C.    SYSTEM  REQUIREMENTS 

System  type 

Memory 

Math  Co-processor  (8087) 

Graphics  card 


Graphics  program 


IBM    PC    or    compatible 
2  5  6K ,     minimum 
Recommended,  not 

required . 

Recommended,  not 

required . 

X,        Y       graphing       program 
like  GRAPHER,  is 

recommended,  not 

required . 

D.    CONCEPT  DESCRIPTION 

1 .   The  Horse  Shoe  Vortex 

This  program,  which  models  the  flow  by  a  series  of 
horse  shoe  vortices,  distributed  over  the  wing,  is  an 
adaptation  of  the  VORLAT  program  by  Moran  [Ref .1] .  He 
develops  a  program  to  find  the  strengths  of  a  series  of 
horse  shoe  vortices  associated  with  a  flat  rectangular  wing. 
The  VORLAT  program  has  its  foundation  in  two-dimensional 
airfoil  theory,  and  Moran  [Ref.l]  adapted  the  theory  to 
wings  of  finite  aspect  ratio. 

The  Users  Manual  will  present  a  short  description  of 
the  VORTEX  program.  For  coverage  of  the  VORLAT  program,  the 
foundation  of  the  VORTEX  program,  consult  Moran  [Ref.l]. 
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a.   Down-wash  and  the  Relationship  to  Flow  Tangency 

The  velocity  induced  at  a  point  in  the  xy  plane 
by  a  horse  shoe  vortex,  also  located  in  the  xy  plane,  can  be 
determined  from: 


Ar 

w  (x,y)  =   

V     J         4n(y-y  ) 


1  + 


V(x-x  f+iy-y  Y 


(x-x  ) 
a 


AT 


4n(y-y  ) 


1  + 


V(x-x/+(y-yf 


(x-x  ) 

a 


7.1 


if  (x    =  x~) , 


AT  / 

w   Uj)  =  — 


4n   \y-yQ      y-yb 


1  .  2 


w    is  the  velocity  induced  at  (x,y),     called  a  control  point. 
AT  is  the  strength  of  the  horse  shoe  vortex,  and  (xa,ya)     and 
(xa>yb)      are   t^le   corners   of   the   horse   shoe   vortex.     The 
center  of  the  horse  shoe  vortex  is  a  field  point. 

The  velocity  at  a  point  on  a  flat  wing,  which 
is  at  an  angle  of  attack  a,  has  components  Vo0cos(a)  and 
VgoS  in  (a)     as  indicated  in  Figure  7-1. 


Vcy  slnCd? 


V^  cos  Ceo 


FLAT    VING 


Velocity  Components  on  a  Flat  Wing 
Figure  7-1 
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In  order  to  ensure  the  flow  is  tangent  to  the 
wing,  the  induced  velocity  or  down-wash-1-  at  that  point  on 
the  wing  must  also  be  Vaosin(a).  Each  vortex  on  the  wing 

contributes  to  the  down-wash^  at  every  control  point  on  the 
wing.  So,  an  equation  for  each  control  point,  as  a  function 
of  all  the  field  points  can  be  written.  An  example  of  the 
equation  for  the  control  point  (Xp,yp)     is. 


-  Y  Ar 


>  Ar  w   (x  ,y  )  =  0 
—   y  y  p'p 
y 


7  .  3 


When  the  control  points  are  combined,  the  result  is  a  set  of 
N  equations  in  N  unknowns. 

b.   Placement  of  the  Horse  Shoe  Vortex 

The  bound  portion  of  the  vortex  (that  portion 
perpendicular  to  the  onset  flow)  is  placed  at  the  1/4  chord 
point  of  each  grid  element  on  the  wing.  Flow  tangency  is 
evaluated  at  the  3/4  chord  point  of  each  grid  element. 


iDown-wash  is  considered  positive  when  its  direction  is 
downward . 

^The  contribution  will  be  positive  when  the  induced 
velocity  is  downward,  or  negative  when  the  induced  velocity 
i  s  upward  . 
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c.   Satisfaction  of  the  Kutta  Condition 

To  satisfy  the  Kutta  condition,  the  vortex 
strength  must  go  to  zero  at  the  trailing  edge.  For  wings 
and  airfoils  of  zero  thickness,  flow  tangency  enforced  at 
the  trailing  edge  appears  to  be  a  corollary  of  the  Kutta 
condi  t ion . 

2 .   Developing  the  Wing 

a.   Planform  Geometry 

It  is  necessary  to  develop  wing  geometry  from 
aspect  ratio,  taper  ratio  and  sweep  angle.  Aspect  ratio  is 
the  span  divided  by  the  average  chord.  Taper  ratio  is  the 
tip  chord  divided  by  the  root  chord,  and  sweep  angle  is  the 
angle  between  the  leading  edge  and  the  y  axis.  Moran 
[Ref.l]  uses  a  fixed  root  chord  and  varies  the  span  as 
necessary  to  achieve  the  required  aspect  ratio  and  taper 
ratio.  The  VORTEX  program  uses  this  approach,  which  works 
well  . 
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Remote    Velocity 


> 


ct 


b/2 


cr   =   1.0 


Planform  and  Variables  used  in  VORTEX  Program 

Figure  7-2 

Aspect   ratio,   AR ,   is  (b^/S).  Span  (b)       is   the 

distance  from  wing  tip  to  wing  tip,  and  area  (S)  is  the 
total  planform  area  of  the  wing. 

Taper  ratio,  A,  is  the  ratio  of  tip  chord  to  root 
chord . 


A  = 


tip  chord         ct 
root  chord        c 


7.4 


Sweep  angle  delta,  A,  is  the  angle  that  the  leading 
edge  makes  with  a  line  drawn  perpendicular  to  the  root 
chord . 

b.   Grid  Geometry 

Using  aspect  ratio,  taper  ratio,  and  sweep 
angle,  the  subroutine  SET85  develops  a  grid.  The  program 
stores  the  coordinates  of  all  necessary  points  in  the 
matrices  WING  and  SECTN.  Results  of  the  program  are  also 
stored  in  WING  and  SECTN  which  are  written  to  data  files 
when  the  program  finishes . 
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In  a  manner  analogous  to  that  used  by  Hough  [Ref.3], 
the  tip  vortices  are  inset  to  improve  accuracy  of  the 
results.   The  inset  distance  is  one  element  wide. 

The  model  keeps  matrix  size  and  the  number  of  grid 
elements  small  to  keep  the  time  necessary  for  solution 
within  reason.  The  model  concentrates  grid  elements  in 
areas  where  the  rate  of  change  is  rapid,  or  where  the  values 
are  most  important.  Those  areas  are  along  the  wing 
boundaries.  An  excellent  method  for  concentrating  the 
points,  while  keeping  the  total  number  of  points  at  a 
minimum,  is  through  cosine  spacing.  Figure  7-3  contains  a 
sample  of  the  layout.  Notice  that  the  rectangular  elements 
do  not  model  the  planform  exactly,  but  in  the  limit  the 
difference  will  be  negligible. 
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Grid  Model  with  Cosine  Spacing 
Figure  7-3 
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3 .   Matrices  Used  in  VORTEX  Program 

The  program  uses  two  identical  matrices  of 
influence  coefficients,  [A]  and  [ AA ] .  [A]  is  used  to  solve 
for  the  AT  of  the  flat  wing,  and  [  AA  ]  is  used  to  solve  for 

o 

the  wing  shape  of  the  elliptically  loaded  wing. 


'A, 
-r  )  =  {AD       (Flat  Wing) 


dx 


rn 


7  .5 


AA 


.Ml  = 


dz 
dx 


(Elliptically  Loaded  Wing) 


7.6 


The  SECTN  matrix  contains  values  for  the  wing 
sections,  (chord,  and  width)  as  well  as  some  of  the  final 
results.  Each  row  of  the  SECTN  matrix  corresponds  to  a 
wing  section  in  the  semi-span.  The  Output  Variables 
section  (H.l.a)  describes  the  columns. 

The  WING  matrix  contains  wing  coordinates  required 
by  the  program  as  well  as  some  of  the  final  results.  Each 
row  of  the  matrix  corresponds  to  an  element  in  the  wing 
semi-span.  The  Output  Variables  section  (H.l.b)  describes 
the  c  o lumns . 


^The   method   used   to   develop   wing   shape   will   be 
discussed  in  section  VII. E. 
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4 .  Solution  of  the  Set  of  Equations 

The  subroutine  GAUSS  solves  the  N  equations  in  N 
unknowns.  The  subroutine  destroys  [A]  in  the  process  but 
returns  the  solution  vector  (AT)  in  place  of  the  right  hand 
sides . 

5 .  Finding  the  Forces  on  the  Elements 

The  Kut ta - Joukowski  theorem  states  that  the  force 
per  unit  span  acting  on  an  element  is. 


—  =  pV~*x  AT 

Ay       «ff 


7.7 


In  this  equation,  ^eff  ^s  the  local  effective  velocity  at 
the  center  of  the  element.  ^eff  ^s  defined  in  equation  7.8 
and  Figure  7-4.  p  is  the  density  and  AT  is  the  incremental 
circulation  around  the  element.  This  circulation  (AT)  is 
the  same  AT  contained  in  equation  7.1.  AT  can  be  termed  the 
incremental  circulation  that  occurs  over  the  element. 
Throughout  this  section,  certain  approximations  and 
dimensional  simplifications  will  be  made.  The  first  is  the 
small  angle  approximation,  where  the  sine  is  approximately 
equal  to  the  angle  in  radians,  and  the  cosine  is 
approximately  equal  to  1.0.  This  approximation  requires 
that  angle  of  attack  for  the  flat  wing  be  small,  generally 
less  than  10  degrees.  The  approximation  also  requires  that 
any  wing  slope  on  the  elliptically  loaded  wing  be  small. 
This   requirement   is   met   by   restricting   the   desired   lift 
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coefficient  to  values  less  than  about  0.5.  In  the  process 
of  determining  coefficients  of  lift,  drag  and  moment, 
certain  dimensional  reference  quantities  arise.  These 
quantities  are  density  (p),  remote  velocity  (V^),  planform 
area  (S),  and  average  chord  (c).  In  performing  dimensional 
analysis  an  arbitrarily  unit  value  can  be  assigned  to  these 
t e  rms . 

a.   Finding  ACp 

The  distribution  of  ACp  over  the  wing  is 
desired,  so  the  force  on  each  element  must  be  converted  to  a 
dimens ionles s  pressure  difference  coefficient.  V  eff  ^s 
de  fined  as . 


>  *  (    r  \  $ 

V  ,,    =  V  cos  a  i    +      V  sin  a  —  w     j 
eff  o°  \     *  ) 


1  .8 


Veff 


Components  of  Effective  Velocity  (Veff) 
Figure  7-4 
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It  is  also  significant  to  note  that  AT  can  be  written  as, 


AF  =  Ar  j 


7.9 


since  the  bound  portion  of  the  horse  shoe  vortex  is  aligned 
with  the  y  axis.  The  result  of  the  force  cross  product, 
Equat  ion  7.7,  is . 


Ay 


w  —  V   sin  a  )AP  i    +  V  cos  a  A  T  k 


7.10 


Setting  p     =  1,  Vco     =     1,   and  incorporating  the  small  angle 
approximation. 


A^ 
Ay 


w  -  a  Ar;  +  &rk 


i .  li 


This  force  is  related  to  ACp .  ACp  is  a  scalar,  rather  than 
a  vector  quantity,  so  the  program  uses  the  force  component 
normal  to  the  wing  to  compute  ACp .   That  component  is. 


AF  =  AjyAf 


7.  12 


Making        the        force        non-dimensional        gives        a        pressure 
difference    coefficient,     ACD. 


AC  =   - 

p       1 


AF 


-pV  AxAv 
2       °° 


7.13 
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Finally,  after  setting  p    -  1  and  V^  =  1,  ACp  for  an  element 
is  . 


2&F        2Ar 

AC  =  =  

P       AjAj    A* 


7.14 


This  provides  the  ACp  for  each  element,  which  can  be  plotted 
versus  chord,  for  each  span-wise  section. 

b.   Finding  the  Lift,  Drag  and  Moment 

Lift,  drag  and  moment  are  also  found  from  the 
vector  force  on  each  element.  Equation  7.10  showed  that  the 
force  on  an  element  has  components  normal  and  tangent  to  the 
wing.  Once  more  setting  p  —  1,  and  Vro  =  1,  the  vector 
force  on  an  element  is. 


W  = 


w  —  sina  lAFAjyi    +  cosaATAv^ 


7  .15 


These  forces  are  oriented  relative  to  the  wing  coordinate 
system,  where  i  is  parallel  to  the  wing  and  k  is 
perpendicular  to  the  wing.  Lift  and  drag  are  normal  and 
tangent,   respectively,   to   the   remote  velocity  (V«>)  ■  The 

wing  is  at  an  angle  of  attack  (a)  to  that  remote  velocity. 
Using  this  information,  the  program  transforms  the  forces  to 
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a  coordinate  system  oriented  to  the  remote  velocity.  For  an 
individual  element  the  lift  and  drag  forces  are; 

AL  =  Elemental  Lift  Force  —  (cos  a)  (Ay AD  (cos  a)  —  (w  —  sin  a)  (AjAD  (sin  a) 
AD  =  Elemental  Drag  Force  =  (cos a MA^AT)  (sin  a)  +   (w  —  sin  a)  (AyAD(  cos  a) 

using  a  small  angle  approximation,  7.16 

AL  -  Elemental  Lift  Force  =  Ay  AT  —  i^aA^AT  =  AjAr 
AD  =  Elemental  Drag  Force  =  AyATa  +  (w  —  a)AyAV  -   wAyAr 

The  total  lift  and  total  drag  acting  on  the  wing  is  a 
summation  of  the  elemental  lifts  and  drags.  Using  wing 
symme  try ; 


L  =  Total  Lift  =  2  AL 

semi  -  span 

7.17 
D  =  TotalDrag  =  2       ^        AD 

semi  —  span 


The  moment  about  the  y  axis   for  an  individual 
element,  with  the  accepted  sign  convention,  is; 

AM  =  Elemental  Moment  =  —  (normal  wing  force)  (x  position  of  center  of  element)  -j     ^3 


AM  =  Elemental  Moment  =  ( -  AyAD  (x      ) 

cen 
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Total  moment  about  the  y  axis  is  a  summation  of 
the  elemental  moments.   Once  again,  using  wing  symmetry. 


M  -  Total  Moment  =  2        ^_        AM 
semi  —span 


7.19 


6 .   Finding  the  Coefficients 

Non - dimens ionaliz ing  the  total  lift  and  drag  gives, 


C,  = 


L         1    2 

2   ■ 


7.20 


C 


D 


o       1 


-p\TS 
2   " 

where  S  -  span  x  average  chord  =  (b  x  c)  ,  p    =  1  and  Vco  =  1. 

Non- dimens ional iz ing   the   total   moment   about   the  y 

axis  gives , 


Mo         ] 


M 


1  .21 


pV~cS 


where  c  is  the  average  chord. 

The  aerodynamic  center  Xac  is  important.  It  is  the 
chord-wise  position  about  which  the  flat  wing  has  zero 
moment.   The  program  finds  Xac  from. 


—  total  moment 

X     =  

K  total  lift 


1  .22 
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This  assumes  that  total  moment  is  taken  about  the  y  axis,  or 
x=0  .  The  program  finds  the  moment  on  the  elliptically 
loaded  wing,  and  uses  the  Xac  from  the  flat  wing  as  the 
reference  axis  to  find  Cj^ac  . 

E.   TECHNIQUE  FOR  SPECIFIED  LOADING 

1 .   Finding  Camber  and  Twist  from  Loading 

Elliptic  loading  in  a  span-wise  direction  results  in 
minimum  induced  drag,  and  loading  is  proportional  to 
circulation.  At  sufficiently  high  aspect  ratio,  a  flat  wing 
with  elliptic  area  distribution  generates  elliptic  loading 
but  weighs  more  than  a  straight  wing.  Wings  with  straight 
leading  and  trailing  edges  also  cost  less  to  manufacture 
than  elliptic  wings.  So,  a  span-wise  elliptic  lift 
distribution  is  generated  by  slightly  twisting  the  wing. 
The  direct  relationship  between  wing  shape  (camber  and 
twist)  and  circulation  distribution  makes  this  generation 
possible  . 
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canber 


twist 


Camber  and  Twist  on  the  Elliptically  Loaded  Wing 

Figure  7-5 

For  a  thin  wing,  elliptic  load  distribution  in  the 

chord-wise   direction   requires   an   approximately   parabolic 

camber  shape.    The  program  specifies  an  elliptic  chord-wise 

load  distribution  and  then  finds  the  associated  shape. 

2  .   Specifying  the  Lift  Coefficient 

Specifying  the  form  and  amplitude  of  the  ideal  lift 
distribution  over  the  wing  determines  the  form  and  amplitude 
of  the  wing  shape.  A  load  distribution  of  elliptic  form  is 
imposed  in  this  case  and  the  corresponding  amplitude  is 
fixed   by   the   desired   ideal   wing   lift   coefficient,   Clj_  . 

3  .   Maximum  Desired  Lift  Coefficient 

The  small  angle  approximation  breaks  down  when  large 
lift  coefficients  are  specified.  To  prevent  this 
occurrence,  the  maximum  desired  lift  coefficient  is 
restricted  to  0.5. 
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4 .  Achieving  the  Desired  Lift  Coefficient 

The  program  scales  a  reference  distribution  of  ACp  , 
which  is  elliptic  in  form,  to  achieve  the  desired  total  lift 
coefficient  for  the  wing. 

5 .  Forces.  Moments  and  Coefficients 

Once  the  ACp  distribution  is  scaled  for  the 
elliptically  loaded  wing,  the  program  finds  forces,  moments 
and  coefficients  in  the  same  manner  as  for  the  flat  wing. 
The  only  difference  is  in  the  coordinate  transformation  used 
to  convert  forces  on  the  wing  to  lift  and  drag.  The 
individual  elements  on  the  elliptically  loaded  wing  are  no 
longer  at  a  uniform  angle  to  the  remote  velocity.  This 
change  in  angle  must  be  taken  into  account  when  performing 
the  coordinate  transformation  which  gives  lift  and  drag 
forces  on  individual  elements. 

The  program  finds  moment  coefficient  about  the 
aerodynamic  center,  C^ac ,  from, 


Mac  Mo  Li\ 

C 


1  .  23 


where  X_ac  is  found  from  the  flat  wing,  c  is  the  average 
chord.  C^q  is  the  pitching  moment  of  the  elliptically 
loaded  wing  about  the  y  axis,  and  Cli  is  the  lift 
coefficient  of  the  elliptically  loaded  wing. 
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F.    FLOW  CHART  FOR  MODULES 


VDRTEX 


SET  WING 
PLANFDRM 
CDMPUTE 
COEFFICIENTS 


SET  GRAPHIC 
DISPLAY 


VIEW  GRAPHIC 
DISPLAY 


END 
PRDGRAM 


GRAP85 


GRAPHERJ 


VLAT85   CDAT85    DELAY    SDIG      RNDUP    PLDT 


CDIT85    CDAT85    LDAD85   DELAY     SHAP85    VIEW 


SET85    DNWS85   GAUSS    MULT2    MULT1     MULT3 
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INPUT  VARIABLES 

There  are  eight  input  variables.   They  are: 

1.  Aspect  Ratio.  AR  ,  defined  as  b^/S,  where  b  is  the 
total  span  and  S  is  the  total  area,  can  have  any 
positive  value  less  than  20.0.  Higher  aspect  ratios 
should  use  lifting  line  theory. 

2.  Taper  Ratio.  LAMDA ,  defined  as  Cc/Cr,  where  Ct  is  the 
wing  tip  chord  and  Cr  is  the  wing  root  chord,  can  have 
any  positive  value  between  0  and  1. 

3.  Leading  Edge  Sweep.  DELTA,  defined  as  the  angle  in 
degrees  that  the  leading  edge  of  the  wing  makes  with  a 
line  perpendicular  to  the  remote  velocity,  can  have 
any  positive  value  between  0  and  60. 

4.  Angle  of  Attack.  ALPHA,  defined  as  the  angle  in 
degrees  made  between  the  flat  wing  and  the  remote 
velocity,  can  have  any  positive  value  between  0  and 
10. 

5.  Number  of  elements  in  a  section.  NX,  defined  as  the 
integer  number  of  elements  in  a  chord  of  the  wing,  can 
have  any  value  between  1  and  10. 

6.  Number  of  elements  in  a  semi-span.  NY,  defined  as  the 
integer  number  of  elements  in  a  semi-span  of  the  wing, 
can  have  any  value  between  1  and  10. 

7.  Desired  lift  coefficient.  CLDSRD,  the  desired  lift 
coefficient  for  the  elliptically  loaded  wing,  can  have 
any  value  between  0  and  0.5. 
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H.   OUTPUT  VARIABLES 
1 .   Tabular  Data 

When  complete,  the  program  writes  four  tabular  data 
files.  One  file,  the  SECTION. MAT  file,  contains  data  for 
each  section  in  the  wing  semi-span.  Another,  the  WING. MAT 
file,  contains  data  for  the  individual  grid  elements  on  the 
wing  semi-span.  The  third,  the  FLAT . DAT  file,  shows  the 
coefficients  for  the  flat  wing.  The  fourth,  the  CAMBER.DAT 
file,  shows  the  coefficients  for  the  cambered  or 
elliptically  loaded  wing. 

a.   SECTION. MAT 

Each  row,  or  line,  in  the  file  represents  a 
section  of  the  wing  semi-span.  Each  column  represents  a 
different  variable  within  the  section. 


Reno~te    Velocity 


> 


chord 


I 


delta    y 


A    Single    section    in    the    Wing    Semi-Span 
Figure    7-6 
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The  columns  are: 

1.  Section  chord  length 

2.  Section  lift  coefficient  per  unit  span  for  the  flat 
wing.  The  integral  of  these  values  from  tip  to  tip  is 
the  total  lift  coefficient. 

3.  Section  drag  coefficient  per  unit  span  for  the  flat 
wing.  The  integral  of  these  values  from  tip  to  tip  is 
the  total  drag  coefficient. 

4.  Section  moment  coefficient  per  unit  span  about  the 
aerodynamic  center  for  the  flat  wing.  The  integral 
of  these  values  from  tip  to  tip  is  the  total  moment 
coefficient  about  the  aerodynamic  center. 

5.  Section  xac  ,  aerodynamic  center  relative  to  x=0  ,  for 
the  flat  wing. 

6.  Section  delta  y 

7.  Section  lift  coefficient  per  unit  span  for  the 
cambered  wing.  The  integral  of  these  values  from  tip 
to  tip  is  the  total  lift  coefficient. 

8.  Section  drag  coefficient  per  unit  span  for  the 
cambered  wing.  The  integral  of  these  values  from  tip 
to  tip  is  the  total  drag  coefficient. 

9.  Section  moment  coefficient  per  unit  span  about  the 
aerodynamic  center  for  the  cambered,  or  elliptically 
loaded  wing.  The  integral  of  these  values  from  tip 
to  tip  is  the  total  moment  coefficient  about  the 
aerodynamic  center. 
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10.  Section   twist   in   degrees   of   the   cambered   wing. 
Positive  angle  of  twist  is  nose  up. 
The  program  arranges  the  rows  from  root  to  tip,  with  the 
root  section  on  the  first  line  and  the  tip  section  on  the 
las  t  line . 

b.   WING. MAT 

Each  row,  or  line  in  the  file,  represents  an 
individual  element  in  the  wing  grid.  Each  column  represents 
a  different  variable  for  that  element.  Figure  7-7  shows  the 
location  of  some  variables  on  a  representative  grid  element. 


Remote  Velocity 


(x,yb) 

(xp,yp)  (force) 
(xp,yp)  (shape) 

(x,ya) 


delta  x 


A  Single  Element  in  the  Wing  Semi-Span 
Figure  7-7 
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The  columns  are: 

1.  Sequential   integer   identifier   of   the   element,   also 
called  IJ  in  the  program. 

2.  Integer  y  position  of  the  element.    1  is  at  the  wing 
root,  NY  is  at  the  tip. 

3.  Integer   x   position   of   the   element.     1   is   at   the 
leading  edge,  NX  is  at  the  trailing  edge. 

4.  x  position  of  center  of  element. 

5.  y  position  of  center  of  element. 

6.  x  position  of  horse  shoe  vortex,  a  field  point. 

7.  Blank,  or  zero. 

8.  ya  position  of  one  corner  of  horse  shoe  vortex. 

9.  yb  position  of  one  corner  of  horse  shoe  vortex. 

10.  xp  trailing  edge  of  the  element,  a  control  point  used 
in  computing  wing  shape  or  AT. 

11.  xp   center   of   the   element,   a   control   point  used   in 
computing  the  forces  acting  on  the  wing. 

12.  yp  center  of  the  element,  a  control  point. 

13.  Delta  x,  chord-wise  dimension  of  individual  element. 

14.  Delta  y,  span-wise  dimension  of  individual  element. 

15.  Flat  wing  incremental  circulation  (AT). 

16 .  Flat  wing  ACp . 

17.  Cambered  wing  incremental  circulation   (AT),   adjusted 
for  C^ef . 

18 .  Cambered  wing  ACp 

19.  Cambered  wing  slope   (3z/9x). 


20.  Cambered  wing  height  (z). 

21.  Cambered  wing  incremental  circulation  (AT),   if  Cj^ref 
were  =  1.0. 

22.  Down-wash  (w)  at  center  of  flat  wing  element  due  to 
trailing  filaments. 

23.  Down-wash  (w)  at  center  of  cambered  wing  element  due 
to  trailing  filaments. 

Some  columns  are  identical.  Separate  columns  were  used  for 
each  variable  during  the  development  phase  to  promote  ease 
in  modification.  The  rows,  or  lines,  of  the  WING. MAT  file 
are  arranged  as  follows. 


Remote  Velocity 


> 


/ 

\ 
\ 

h 

2 

3        \ 

Grid  Element  Numbering  on  the  Semi-Span 
Figure  7-8 

The   first   row   is   the   leading   edge   element   in   the   root 

section.   The  next  element  is  immediately  behind  the  leading 

edge  and  so  on.  The  last  row  is  the  trailing  edge  element  in 

the  wing-tip  section. 
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c.   FLAT. DAT 

This  file  contains  coefficients  and  results  for 
the  flat  wing.   An  example  is  shown. 


FLAT  WING 

CL 

= 

.177721 

CD 

= 

.006516 

CD/CL2 

= 

.206312 

CMAC 

= 

.000000 

XAC 

= 

.279935 

AR 

= 

1 

.330000 

LAMDA 

= 

.500000 

DELTA 

= 

25 

.000000 

NX 

= 

8 

NY 

= 

8 

ALPHA 

= 

5 

.000000 

CLDSRD 

= 

.200000 

ELLIP      =  YES 

CLa,    Lift   Curve   Slope 

per/Deg=  .035544 

d.        CAMBER.DAT 

This  file  contains  coefficients  and  results  for 
the  cambered,  or  elliptically  loaded,  wing.  An  example  is 
shown . 

CAMBERED  WING 


CL 

= 

.199311 

CD 

= 

.009040 

CD/CL2 

= 

.227557 

CMAC 

= 

- 

.062732 

AR 

= 

1 

.330000 

L^MDA 

= 

.500000 

DELTA 

= 

25 

.000000 

NX 

= 

8 

NY 

= 

8 

ALPHA 

= 

5 

.000000 

CLDSRD 

= 

.200000 

ELLIP 

= 

YES 
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I .   SAMPLE  PROBLEM 

A  sample  problem  will  illustrate  use  of  the  VORTEX 
program.  A  wing  planform  with  the  following  characteristics 
will  be  analyzed. 


Aspec  t  ratio 
Taper  ratio 
Leading  Edge  Sweep 
Angle  of  Attack 
Number  of  span  elements    8 
Number  of  chord  elements   8 
Desired  Lift  Coefficient   0.2 
The  planform  looks  like  Figure  7-9 


1.  333 

0.  5 

2  5  degrees 

5  degrees 


f—  0.5000  — | 


Reno~te    Velocity 


> 


0.5000 


1.0000 


Planform  used  in  the  Sample  Problem 
Figure  7-9 
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1.   Starting  the  Program 

After  turning  the  computer  on  go  to  the  DOS  prompt, 
which  generally  looks  like  this. 
C:> 

Change  the  program  to  the  VORTEX  directory  by  typing 
CD\VORTEX  [return] 

the  screen  should  now  look  like. 
C\VORTEX:> 

To  start  the  program,  type 
VORTEX  [return] 

The   program   will   start   and   the   following   menu   will   be 
displayed . 


MAIN  PROGRAM  MENU 

1.  Set  Wing  Planform,  compute  Coefficients  and  Shape 

2.  Set  Graphic  Display 

3.  View  Graphic  Display 

4 .  End  the  Program 


Use  this  option  to 

set  up 

yoTir 

wing'  plant orm. 

It 

will 

also  compute  all  the 

coef  f 

icients  you  need , 

Lift, 

Drag,  Moment, 

as 

well 

as  the  shape , 

if  you 

want 

an  elliptic  load 

disti 

-ibution  . 

I 

Main  Program  Menu 
Figure  7-10 

2 .   Main  Program  Menu 

There  are  four  options  in  the  MAIN  PROGRAM  MENU. 
All  coefficients  and  results  are  computed  in  option  1. 
Whenever  the  planform  is  changed,  option  1  must  be  selected 
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to  recompute  the  coefficients.  The  coefficients  and  results 
are  saved  to  data  files  after  being  computed  in  option  1. 
The  data  files  allow  you  to  use  options  2  and  3  as  often  as 
desired  without  re-computing  the  results  for  the  planform. 

Help  menus  are  available.  Press  Fl  to  toggle  the 
help  menus  on  or  off. 

3 .   Planform  and  Coefficient  Menu 

The  program  saves  planform  variables  and  results 
from  the  most  recent  solution  and  shows  them  in  the  PLANFORM 
and  COEFFICIENT  MENU.   The  menu  looks  like  Figure  7-11. 


PLANFORM  AND  COEFFICIENT  MENU 


1. 
2. 

3  . 

4  . 
5. 
6. 

7  . 

8  . 

9  . 
10. 


Aspect  Ratio  , 

Taper  Ratio  , 

Leading  Edge  Sweep  Angle  (degrees)  .  .  , 

Ancle  of  Attack  (degrees)  , 

Number  of  Elements  in  a  Section  , 

Number  of  Elements  in  a  Semi -Span  . . . 

Compute  Wing  Shape  with  Load 

E'esir-ed  Lift  Coefficient  , 

Compute  Cofficients  for  this  Planform 
Return  to  Main  Menu 


1 . 330000 
. 500000 
25.000000 
5.000000 
8 
8 
YES 

.200000 


Aspect  ratio  has  values 
between  0  and  10. 


Planform  and  Coefficient  Menu 
Figure  7-11 

Options  1,  2,  and  3  change  Aspect  Ratio,  Taper 
Ratio,  and  Leading  edge  sweep  angle,  respectively. 

Option  4  changes  the  angle  of  attack. 

Option  5  changes  the  number  of  chord-wise  grid 
elements ,  and 
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Option  6  changes  the  number  of  span-wise  grid 
elements  in  a  semi-span. 

Option  7,  Compute  Wing  Shape  with  Load,  is  YES  or 
NO.  Because  the  time  required  to  compute  wing  shape  is 
extensive,  that  option  may  not  be  desirable  for  all 
prob lems . 

Option  8  changes  the  desired  lift  coefficient,  C^ 
of  the  elliptically  loaded  wing. 

Option  9  computes  the  coefficients.  The 
coefficients  are  written  to  four  files,  SECTION. MAT, 
WING. MAT,  FLAT. DAT,  and  CAMBER.DAT.  Section  H  of  this 
thesis  contains  a  description  of  the  elements  in  the  files. 

Option  10  returns  to  the  Main  Program  Menu  to  view 
the  re  suits . 
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4 .   Graphics  Program  Menu 

Choosing  selection  2  of  the  Main  Program  Menu 
displays  the  Graphics  Program  Menu.  An  example  is  shown  in 
Figure  7-12. 


GRAPHICS  PROGRAM  MENU 

Type  of  Wing  

Section  Number  

Return  to  Main  Menu 


Spanwise  Plots 

4.  d(CL)/dy 

5.  d(CD)/dy 

6.  d(CMAC)/dy 

7 .  Camhex-ed  Wing  Twi.s 

Chordwise  Plot 

8.  Delta  Cp 

9 .  Wing  Shape 


CAMBERED 
5 


Use  th 

is 

option 

to 

return 

to  the 

Ma 

in  Pro 

gram  Mem:, 

where 

you 

can  view 

the 

graph 

you 

have 

set 

up 

here . 

Graphics  Program  Menu 
Figure  7-12 

Option  1  changes  the  type  of  wing  that  will  be 
displayed,  either  flat  or  cambered. 

Option  2  changes  the  section  of  the  wing  which  is 
plotted  with  options  8  and  9.  . 

After  setting  options  1  and  2,  select  the  desired 
plot  from  options  4  through  9.  The  necessary  data  files 
will  be  prepared  and  the  program  will  return  to  the  Main 
Program  Menu,  where  the  selected  plot  can  be  seen  using  the 
View  Graphic  Display  option.   Examples  follow. 
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CL 

FLAT   WING 
=      .1693,   AOA   =      5.00    deg 

0.30  - 

0.24  - 

£0.18- 

o 

^o  0.12  - 

0.06  - 

0.00  -\ —    — 
0.0 

i                  i                  i                   i 
0.1                0.2               0.3               0.4 

SPAN-WISE   COORDINATE   Y 

0.5 

FLAT  WING  d(CL)/dy  vs  SPAN 
Figure  7-13 


CAMBERED   WING 

Cli      =      .1993,    CLi    =         .20 


0.30  -i 


0.24  - 


£0.18  H 


^o0.12  J 


0.06  - 


0.00 


0.0  0.1  0.2  0.3  0.4 

SPAN-WISE    COORDINATE   Y 


0.5 


CAMBERED  WING  d(CL)/dy  vs  SPAN 
Figure  7  - 14 
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FLAT   WING 

CD 

=      .0062,   AOA   =      5.00    deg 

0.020  - 

0.012  - 

$    0.004  - 

\ 
/ — - 

Q 

O                  0. 

0 

0.1               0.2              0.3              0.4 

0.5 

-a  -0.004  - 

-0.012  - 

-0.020  - 

SPAN-WISE    COORDINATE   Y 

FLAT  WING  d(CD)/dy  vs  SPAN 
Figure  7-15 


CAMBERED    WING 

CD         =      .0090,    CLi    =         .20 


0.020  -i 


0.016  - 


-^0.012  - 


O 
o 

^u  0.008 


0.004  H 


0.000 


0.0  0.1  0.2  0.3  0.4 

SPAN-WISE    COORDINATE  Y 


0.5 


CAMBERED  WING  d(CD)/dy  vs  SPAN 
Figure  7-16 
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FLAT   WING 
CMAC    =      .0000,   AOA   =      5.00    deg 

0.0110  - 

0.0070  - 

\ 

>N 

\^ 

\    0.0030  - 

< 
2 

\ 

•B-o.ooioOr 

0              0.1               0.2  \        0.3              0.4 

i 
0.5 

-0.0050  - 

-C.0090  - 

SPAN-WISE   COORDINATE  Y 

FLAT    WING    d(CMAC)/dy    vs     SPAN 
Figure    7-17 


CAMBERED   WING 

CMAC    =    -.0625.    CLi    =         .20 


-0.000 


0. 


-0.018  - 


-0.036  - 


-0.054  - 


-0.072  - 


-0.090  J 


0.1  0.2  0.3  0.4 


0.5 


SPAN-WISE    COORDINATE   Y 


CAMBERED  WING  d(CMAC)/dy  vs  SPAN 
Figure  7-18 
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2.0  - 
1.6  - 

&1.2- 

_i 

UJ 

Q0.8  - 
0.4  - 
n  n  - 

FLAT  WING,    SECTION         1 

LE                                                                                                                     IS. 

U.U    1                           i                            |                            |                            |                           | 

0.0              0.2               0.4              0.6               0.8               1.0 

CHORD-WISE    COORDINATE    X 

FLAT  WING  DELTA  CP  vs  CHORD,  Section  1 
Figure  7-19 


CAMBERED   WING,    SECTION         1 


0.30  -i 


o0.18 


0.0  0.2  0.4  0.6  0.8 

CHORD-WISE    COORDINATE    X 


CAMBERED  WING  DELTA  CP  vs  CHORD,  Section  1 

Figure  7-20 
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FLAT   WING,    SECTION         4 


2.0  - 

LE. 

JS. 

1.6  - 

O   '-2  - 

P 

DELTA 

bo 

i 

0.4  - 

0.0  - 

I                        1 

I                     T*  '                  ! 

0.0  0.2  0.4  0.6  0.8 

CHORD-WISE   COORDINATE   X 


1.0 


FLAT  WING  DELTA  CP  vs  CHORD,  Section  4 
Figure  7-21 


CAMBERED    WING,    SECTION         4 


0.30  - 

J-E. 

J£. 

0.24  - 

Cj0.1S  - 

DELTA 

o 

0.06  - 

0.00  - 

I                            I 

i                          I                          I 

0.0  0.2  0.4  0.6  0.8 

CHORD-WISE   COORDINATE   X 


1.0 


CAMBERED  WING  DELTA  CP  vs  CHORD,  Section  4 

Figure  7-22 
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FLAT  WING.    SECTION         7 

2.0  -i 

X.E.                                                                 IX 

1.6  - 

U  1-2 

|                         | 

3 

UJ 

Q  0.8  - 

0.4  - 

0.0  - 

0. 

H T                     I                   '  1 

0               0.2               0.4               0.6               0.8 

CHORD-WISE    COORDINATE    X 

i 
1.0 

FLAT  WING  DELTA  CP  vs  CHORD,  Section  7 
Figure  7-23 


CAMBERED   WING.    SECTION 


0.30 


0.24  - 


&0.18  - 


UJ 

a  0.12  - 


0.06  -1 


0.00 


0.0  0.2  0.4  0.6  0.8 

CHORD-WISE    COORDINATE    X 


1.0 


CAMBERED  WING  DELTA  CP  vs  CHORD,  Section  7 

Figure  7-24 
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FLAT   WING.    SECTION         8 

2.0  - 

L.E.                                                                  Ji 

1.6  - 

1,,. 
O  1.2  - 

< 

\— 
_J 
UJ 
Q  0.8  - 

(            ' 

0.4  - 

L_ 

o.o  H 

0. 

0               0.2               0.4              0.6               0.8 

CHORD-WISE   COORDINATE   X 

i 
1.0 

FLAT  WING  DELTA  CP  vs  CHORD,  Section 
Figure  7-25 


CAMBERED   WING,    SECTION         8 


0.30  -i 


IX 


I.E. 


0.8 
CHORD-WISE   COORDINATE   X 


1.0 


CAMBERED  WING  DELTA  CP  vs  CHORD,  Section 

Figure  7-26 
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0.030  - 

M 

0.010  - 
Q 

ac 

8 

u  -0.010  - 

h- 

X 

g 

UJ 

I  -0.030  - 
o 

X 

5  -0.050  - 

-0.070  - 

CAMBERED   WING,    SECTION         1 
vertical    scale    exaggerated 

LE.                                                                                                                     J.E. 

Di              0.2          ^0^4^^      °-'6               °-'8               1'-'° 
CHORD-WISE    COORDINATE   X 

CAMBERED  WING  SHAPE  vs  CHORD,  Section  1 
Figure  7-27 


CAMBERED   WING,    SECTION         4 
vertical    scale   exaggerated 
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CAMBERED  WING  SHAPE  vs  CHORD,  Section  4 
Figure  7-28 
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CAMBERED    WING.    SECTION         7 
vertical    scale    exaggerated 
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CAMBERED  WING  SHAPE  vs  CHORD,  Section  7 
Figure  7-29 


CAMBERED   WING,    SECTION         8 
vertical    scale    exaggerated 
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CAMBERED  WING  SHAPE  vs  CHORD,  Section 
Figure  7-30 
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CAMBERED   WING 
CLi      =         .20 
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3.2  - 
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a: 

024- 
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CAMBERED  WING  TWIST  vs  SPAN 
Figure  7-31 
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J.   METHOD  FOR  CHANGING  LOADING  DISTRIBUTION 

This  is  an  optional  section  and  may  be  skipped  without 
loss.  Elliptic  loading  is  generated  in  the  SET85 
subroutine.  The  elements  which  contain  the  specified 
loading  are  WING(IJ,21).  If  you  wish  to  change  the  loading 
from  elliptic,  realize  that  the  grid  elements  are  not 
uniform  in  size. 

When  generating  a  ACp  function,  consider  the  shape.  The 
load  at  the  wing  tips  and  trailing  edge  must  go  to  zero. 
The  load  at  the  leading  edge  should  preferably  be  zero  to 
avoid  a  singularity  at  that  point. 

If  you  change  any  of  the  subroutines,  the  program  must 
be  re-compiled  with  a  FORTRAN  compiler  to  include  the 
changes  into  the  VLAT85.EXE  file.  The  original  version  of 
the  program  was  compiled  by  a  MICRO-SOFT  4.2  compiler  using 
the  following  commands. 
FL/FPc  VLAT85.FOR 
and 

FL/FPc  GRAP85.FOR 

If  you  have  a  math  co  -  proces  sor  ,   the  FPc  option  allows  the 
program  to  use  the  co-processor  but  does  not  require  one. 

K.   METHOD  FOR  CHANGING  THE  NUMBER  OF  GRID  ELEMENTS 

This  is  also  an  optional  section  and  can  be  skipped.  In 
its  original  form,  the  program  works  with  up  to  100  grid 
elements  in  the  wing  semi-span.  The  program  also  has  a 
maximum  of  10  chord  sections.    Fewer  grid  elements  or  chord 
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sections  will  run  without  modification.  If  you  desire  more 
than  100  grid  elements  or  10  chord  sections,  you  must  modify 
the  program.  The  changes  are  not  extensive,  but  require 
that  the  programs  be  recompiled,  as  described  in  section  J 
above . 

Increase  grid  size  in  square  increments.  Modify  the 
matrix  elements  as  follows.  N  is  the  square  dimension  and 
corresponds  to  NX  or  NY. 

For  example,   if  you  want  12  wing  sections   or  12  grid 
elements   in  each  chord  section,   make  the   grid  12x12   (144 
elements).   In  this  example,  N  =  12,  and  (N*N)+1  =  145.   The 
matrices   that   must   be   changed   are   listed   below.     The 
required  dimensions  are  shown. 
A( (N*N)+1 , (N*N)+1) 
AA( (N*N)+1 , (N*N)+1) 
AS ( (N*N)+1 , (N*N)+1) 
WING( (N*N)+1 , 23) 
SECTN(N+1 , 10) 


lOf 


These   are   the   line   numbers   in   the   programs   which 
contain  matrix  variables  that  must  be  changed. 

Var  iable 


Program 

A() 

AA() 

AS() 

WING( ) 

SECTN( 

VLAT85 

76 

77 

77 

74 

74 

LOAD85 

56 

57 

57 

54 

54 

SHAP85 

54 

55 

55 

52 

52 

GAUSS 

19 

- 

- 

- 

- 

MULT1 

19 

20 

20 

17 

17 

MULT2 

19 

20 

20 

17 

17 

MULT  3 

19 

20 

20 

17 

17 

GRAP85 

52 

53 

53 

50 

50 

DNWS85 

- 

- 

- 

23 

23 

SET85 

_ 

_ 

_ 

69 

69 

Before  increasing  the  grid  density,  consider  the 
tradeoff  of  memory  requirement  and  solution  speed.  The  10  X 
10  grid  contains  100  elements  and  100  unknowns.  A  12  x  12 
grid  contains  144  elements  and  144  unknowns.  A  20  x  20  grid 
contains  400  elements  and  400  unknowns.  More  than  4  times 
the  number  of  computations  are  needed  to  solve  a  400  x  400 
versus  a  100  x  100  matrix.  You  should  also  increase 
precision  if  you  increase  the  number  of  grid  elements.  All 
elements  will  be  smaller  and  computer  accuracy  may  be  less 
than  the  distances  between  points  in  the  grid. 
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APPENDIX  A.  SOURCE  CODE  LISTINGS  FOR  VORTEX  PROGRAM 
VLAT8  5 


* 
* 
* 


PROGRAM  VLAT85 

THIS  IS  THE  MAIN  DRIVER  FOR  A  NUMBER  OF  SUBROUTINES  THAT  ARE 
LINKED  TOGETHER  TO  FORM  A  PROGRAM  THAT  COMPUTES  THE  FORCES 
ON  A  PLANFORM  IN  INVICID,  INCOMPRESSIBLE  FLOW. 
THE  SUBROUTINES  CALLED  ARE  LISTED  HERE  AND  AT  THE  END  OF  THE 
PROGRAM  AS  INCLUDE  FILES. 


* 


THIS  PROGRAM  IS  A  HIGHLY  MODIFIED  VERSION  OF  A  PROGRAM  IN  * 

'AN  INTRODUCTION  TO  THEORETICAL  AND  COMPUTATIONAL  AERODYNAMICS'  * 

BY  JACK  MORAN,  WILEY  AND  SONS,  NEW  YORK,  1984,  FOUND  ON  PAGE  * 

151.  * 

* 

PROGRAM  VARIABLES  ARE:  * 

AR      :  ASPECT  RATIO  * 


LAMDA  :  TAPER  RATIO 

ADELTA  :  LEADING  EDGE  SWEEP  ANGLE,  IN  DEGREES 

DELTA  :  LEADING  EDGE  SWEEP  ANGLE,  IN  RADIANS 

NX  :  NUMBER  OF  GRID  ELEMENTS  IN  A  SECTION 

NY  :  NUMBER  OF  GRID  ELEMENTS  ALONG  A  SEMI -SPAN 


* 
* 
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* 
* 
* 

* 

* 

* 

* 

* 
* 

* 
* 


ALPHA   :  ANGLE  OF  ATTACK,  IN  DEGREES 

CLDSRD  :  DESIRED  LIFT  COEFFICIENT  WITH  ELLIPTIC  LOADING 

ELLIP   :  CHARACTER  FLAG  TO  COMPUTE  THE  SHAPE  FROM  AN 

ELLIPTIC  LOADING 
PI     :  IS  PI 

NEQNS   :  NX*NY,  THE  TOTAL  NUMBER  OF  EQUATIONS 

A  (  )   :  A  COEFFICIENT  MATRIX,  USED  FOR  THE  FLOW  TANGENCY  * 

GETS  TRANSFORMED  IN  GAUSS  * 

AA  (  )  :  THE  SAME  AS  A  BUT  NOT  TRANSFORMED  IN  THE  PROGRAM  * 

AS  (  )  :  ANOTHER  COEFFICIENT  MATRIX,  USED  FOR  THE  DOWN -WASH  * 

AT  THE  CENTER  OF  THE  GRID  ELEMENT.   REQUIRED  TO  * 

GET  THE  FORCES  ACTING  ON  EACH  GRID  ELEMENT.  * 

WING(  ):  A  LARGE  ARRAY  HOLDING  VARIOUS  IMPORTANT  COORDINATES  * 

FOR  EACH  GRID  ELEMENT,  AS  WELL  AS  THE  RESULTS  FOR  * 

EACH  GRID  ELEMENT.   A  COMPLETE  DESCRIPTION  OF  EACH  * 

ELEMENT  IN  THE  ARRAY  IS  CONTAINED  IN  THE  SUBROUTINE  * 

SET85  * 

SECTNQ:  A  LARGE  ARRAY  HOLDING  SOME  DIMENSIONS  FOR  THE  * 

SECTIONS,  AS  WELL  AS  THE  COEFFICIENTS  FOR  EACH  * 

SECTION.   A  COMPLETE  DESCRIPTION  OF  EACH  ELEMENT 

IN  THE  ARRAY  IS  CONTAINED  IN  THE  SUBROUTINE  SET85.  * 

UN     :  SET  TO  5,  FOR  KEYBOARD  INPUT  * 

IOUT    :  SET  TO  6  FOR  SCREEN  OUTPUT  * 

JOUT    :  SET  TO  10  FOR  WING  ARRAY  INPUT  AND  OUTPUT  * 

KOUT    :  SET  TO  11  FOR  SECTN  ARRAY  INPUT  AND  OUTPUT  * 

LOUT    :  SET  TO  12  FOR  LAST  PARAMETER  INPUT  AND  OUTPUT  * 
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SET  TO  13  FOR  FLAT  PLATE  WING  COEFFICIENTS 

SET  TO  14  FOR  CAMBERED  WING  COEFFICIENTS 

A  FLAG.   IF  THE  PARAMETERS  ARE  UNCHANGED  FROM  THE 

MOST  RECENT  CALCULATION  OF  COEFFICIENTS,  SET  =  0 

IF  ANY  PARAMETERS  ARE  CHANGED,  SET  =  1 

MENU  INPUT  VARIABLE 

THE  ESCAPE  CHARACTER  (27) 

INTEGER,  1  IF  FLAT,  2  IF  ELLIPTIC 

INTEGER,  SECTION  NUMBER  TO  BE  PLOTTED. 


SUBROUTINES  CALLED  ARE: 

CDAT85  ,  FORTRAN,  GENERAL  DATA  ENTRY,  REAL 

CDIT85  ,  FORTRAN,  GENERAL  DATA  ENTRY,  INTEGER 

LOAD85  ,  FORTRAN,  COMPUTE  LOAD  DISTRIBUTION  FROM  A  GIVEN 

WING  SHAPE. 
SHAP85  ,  FORTRAN,  COMPUTE  SHAPE  FROM  A  GIVEN  LOAD 

DISTRIBUTION. 
DELAY   ,  FORTRAN,  SLOW  DOWN  THE  PROGRAM  SO  ERROR  MESSAGES 

CAN  BE  SEEN  BEFORE  BEING  ERASED. 


* 

MOUT 

* 

NOUT 

* 

IMOD 

* 

* 

* 

ICHOIC 

* 

ESC 

* 

ITYPE 

* 

ISECT 

* 

* 
* 
* 
* 
* 

* 


********************************************************* 
PROGRAM  VLAT85 
REAL  LAMDA 
CHARACTER* 3  ELLIP 
CHARACTERS  ESC 
COMMON  PI,B,WING(101,23) ,  SECTN(11 , 10) ,AR, LAMDA, DELTA, NX.NY, 
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*  ALPHA .CLDSRD.ELLIP.XAC.CLA 
COMMON/COF/A(101 , 101) ,  NEQNS 
COMMON/CAM/AA(101,101) ,  AS (101 ,101) 
COMMON/FILS/IIN , IOUT , JOUT , KOUT , LOUT , MOUT , NOUT 
IMOD  =   1 
AR  =    2.0 
LAMDA  =1.0 
ADELTA=  10.0 
NX  =     4 
NY  =    4 
ALPHA  =5.0 
ELLIP  =  'YES' 
CLDSRD=  0.5 
ITYPE  =  1 
ISECT  =  1 
XAC  =   0.0 
PI  =  2.0*ACOS(0.0) 
ESC  =  CHAR(27) 
********************************************* 

*  * 

*  SET  INPUT/OUTPUT  VARIABLES  AND  OPEN  THE  NECESSARY  FILES         * 

*  MOUT  AND  NOUT  ARE  NOT  OPENED  UNTIL  THE  COEFFICIENTS  HAVE  BEEN   * 

*  COMPUTED .  * 

*  * 
*********************************************************************** 

UN  =  5 

K 


IOUT  =  6 

JOUT  =  10 

KOUT  =11 

LOUT  =12 

MOUT  =13 

NOUT  =  14 

OPEN(UNIT  =  UN) 

OPEN (UNIT  =  IOUT) 

OPEN (UNIT  =  JOUT,  FILE-' WING .MAT' ) 

OPEN (UNIT  =  KOUT,  FILE=' SECTION. MAT' ) 

OPEN (UNIT  =  LOUT,  FILE-' PARAMS . DAT' ) 
10    FORMAT  (I6/3(F6.2,/),2(I3/),F6.2/A/F6.2/I2/I2/F6.2) 
30    FORMAT  (23F10.6) 
40    FORMAT  (10F10.6) 
50     FORMAT  ('  THAT  IS  NOT  A  VALID  OPTION.  ENTER  AN  OPTION  BETWEEN  1' 

*  '  AND  11') 
60    FORMAT  ('  ' ,A,A,F10.6) 
62     FORMAT  ('  ' ,A,A,I3) 
64    FORMAT  ('  ' ,A,A,A) 
************************************************** 

*  * 

*  READ  THE  LAST  SET  OF  PARAMETERS,  IF  UNCHANGED,  READ  WING  AND    * 

*  SECTN  DATA  * 

*  * 
*********************************************************************** 

READ  (LOUT,10,ERR=70,END=70)  IMOD , AR , LAMDA , ADELTA , NX , NY , ALPHA , 
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*  ELLIP,CLDSRD,ITYPE,ISECT,XAC 
70    DELTA  =  ADELTA*PI/180. 
B  =  AR*(1.0+LAMDA)/2.0 
NEQNS  =  NX*NY 
IF  (IMOD.EQ.0)  THEN 

READ  (JOUT.30)  ( (WING(I , J) , J=l , 23) , 1=1 .NEQNS) 
READ  (KOUT.40)  ( (SECTN(I , J) , J-l, 10) , I-l.NY) 
ELSE 
END  IF 

*********************************************************************** 

*  * 

*  CLEAR  THE  SCREEN  AND  WRITE  THE  MAIN  MENU.  * 

*  THIS  USES  THE  ANSI. SYS  CONTROL  CODES  TO  CLEAR  THE  SCREEN  AND  * 

*  POSITION  THE  CURSOR.   IT  ALSO  USES  FLASHUP  WINDOWS  FOR  THE  * 
MENU  PROMPTS.  * 

*  * 
*********************************************** 

100   WRITE  (IOUT,*)  ESC,'[2J' 
WRITE  (IOUT,*)  '"C-ALL/' 
WRITE  (IOUT,*)  '~W=LOAD/' 
WRITE  (IOUT, 60)  ESC , ' [09 ; 69H' , AR 
WRITE  (IOUT, 60)  ESC , ' [ 10 ; 69H' , LAMDA 
WRITE  (IOUT, 60)  ESC , ' [ 11 ; 69H' , ADELTA 
WRITE  (IOUT, 60)  ESC ,'[ 12 ; 69H' .ALPHA 
WRITE  (IOUT, 62)  ESC , ' [ 13 ; 69H' ,NX 
WRITE  (IOUT, 62)  ESC , ' [ 14 ; 69H' ,NY 
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WRITE  (IOUT, 64)  ESC , ' [ 15 ; 69H' , ELLIP 
WRITE  (IOUT, 60)  ESC , ' [ 16 ; 69H' , CLDSRD 
NEQNS  =  NX*NY 
******************************************************* 

*  * 

*  READ  THE  OPTION.   IF  OUT  OF  RANGE,  PRINT  ERROR  MSG ,  AND  START   * 

*  OVER .  * 

*  * 
*********************************************************************** 
110   READ(IIN,*,END=110,ERR=110)  ICHOIC 

WRITE  (IOUT,*)  ESC, ' [ 2J ' 

WRITE  (IOUT,*)  ESC, ' [00;00H' 

IF  ((ICHOIC. LT.l) .OR. (ICHOIC.GT.il))  THEN 

WRITE  (IOUT, 50) 

CALL  DELAY(6) 

GO  TO  100 
ELSE 
ENDIF 
*********************************************************************** 

*  * 

*  MAKE  ALL  CHANGES  TO  PARAMETERS  AS  DESIRED.  * 

*  VALUES  PASSED  TO  CDAT  AND  CDIT  ARE  LOW,  HIGH,  VARIABLE,  AND     * 

*  NAME  OF  THE  VARIABLE.   IMOD  IS  RETURNED  WITH  VALUE  1.  * 

*  * 
*********************************************************************** 

IF  ( ICHOIC. EQ.l)  CALL  CDAT85(0 . 0 , 20 . 0 , AR, 
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*  'ASPECT  RATIO  \IMOD) 
IF  (ICHOIC.EQ.2)  CALL  CDAT85 (0 . 0 , 1 . 0 ,LAMDA, 

*  'TAPER  RATIO  ',IMOD) 
IF  (ICHOIC.EQ.3)  THEN 

CALL  CDAT85(0.0,60.0,ADELTA, 

*  'LEADING  EDGE  SWEEP  ANGLE  IN  DEGREES  \IMOD) 

DELTA  =  ADELTA*PI/180. 
ELSE 
END  IF 
IF  (ICH0IC.EQ.4)  CALL  CDAT85(0 . 0 , 10 . 0 .ALPHA, 

*  'ANGLE  OF  ATTACK  IN  DEGREES  '.IMOD) 
IF  (ICHOIC.EQ.5)  CALL  CDIT85 (0 , 10 ,NX, 

*  'NUMBER  OF  ELEMENTS  IN  A  SECTION      ',IMOD) 
IF  (ICHOIC.EQ.6)  CALL  CDIT85 (0 , 10 ,NY, 

*  'NUMBER  OF  SECTIONS  IN  A  SEMI-SPAN    '.IMOD) 
IF  (ICHOIC.EQ.8)  CALL  CDAT85(0 . 0 , 0 . 5 , CLDSRD , 

*  'DESIRED  LIFT  COEFFICIENT  ',IMOD) 
IF  (ICHOIC.EQ.7)  THEN 

IMOD  =  1 

IF(ELLIP.EQ. 'YES')  THEN 

ELLIP  =  'NO  ' 
ELSE 

ELLIP  =  'YES' 
ENDIF 
ELSE 
ENDIF 
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******************************************************* 

*  * 

*  RUN  THE  MAIN  PART  OF  THE  PROGRAM.   COMPUTE  THE  LOAD,  FROM  THE   * 

*  FLAT  PLATE  SHAPE,  AND  THEN  THE  SHAPE,  FROM  THE  ELLIPTIC  LOAD    * 

*  IF  THE  ELLIP  FLAG  IS  SET.  * 

*  * 
*********************************************************************** 
115    IF  (ICHOIC.EQ.9)  THEN 

IMOD  =  0 

WRITE  (IOUT,*)  ESC,'[2J' 

WRITE  (IOUT,*)  ESC, ' [00;00H' 

WRITE  (IOUT,*)  '~C=ALL/' 

WRITE  (IOUT,*)  '~W=WAIT/' 

CALL  LOAD85 

IF  (ELLIP. EQ. 'YES')  CALL  SHAP85 
ELSE 
ENDIF 
*********************************************************************** 

*  * 

*  OPTION  10,  RETURN  TO  THE  MAIN  MENU.  DO  A  NORMAL  END.  WRITE  THE   * 

*  PARAMETERS  AND  DATA  TO  FILES  AND  CLOSE  THEM.  * 

*  * 
*********************************************************************** 

IF  (ICHOIC.EQ.10)  THEN 
REWIND  (J OUT) 
REWIND  (KOUT) 
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REWIND  (LOUT) 

WRITE (JOUT, 30)  ((WING(I.J) ,J=1,23) ,I=1,NEQNS) 

WRITE (KOUT, 40)  ((SECTN(I.J) ,J-1,10) ,1=1, NY) 

WRITE (LOUT, 10)  IMOD , AR.LAMDA, ADELTA , NX , NY , ALPHA, ELLIP , 
*    CLDSRD,ITYPE,ISECT,XAC 

CLOSE  (JOUT) 

CLOSE  (KOUT) 

CLOSE  (LOUT) 

GOTO  1000 
ELSE 
END  IF 
GO  TO  100 
1000  WRITE  (IOUT,*)  ESC,'[2J' 

WRITE  (IOUT,*)  ESC, ' [00;00H' 
WRITE  (IOUT,*)  '~C=ALL/' 
WRITE  (IOUT,*)  '~W=MAIN/' 
END 

*  * 

*  THESE  ARE  ALL  THE  SUBROUTINES  THAT  ARE  LINKED,  * 

*  * 

$ INCLUDE: 'SET85.FOR' 
$ INCLUDE : ' DNWS  8  5 . FOR ' 
$INCLUDE: 'GAUSS. FOR' 
$ INCLUDE: 'CDAT85.FOR' 


r 


$INCLUDE : ' CDIT85 . FOR' 
$ INCLUDE: 'LOAD85.FOR' 
$ INCLUDE: 'SHAP85.FOR' 
$ INCLUDE: 'DELAY. FOR' 
$ INCLUDE: 'MULT1.FOR' 
$ INCLUDE: 'MULT2.FOR' 
$ INCLUDE: 'MULT3.FOR' 
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LOAD85 


*  * 


* 
* 


* 
* 


SUBROUTINE  LOAD85 


SUBROUTINE  TO  DEVELOP  A  LOADING,  OR  VORTEX  DISTRIBUTION,  THAT  * 

IS  ASSOCIATED  WITH  A  FLAT  PLATE  WING.   THE  LOADING  GENERATED  * 

IS  RELATED  TO  THE  PRESSURE  DIFFERENTIAL,  DELTA  CP,  AND  HAS  THE  * 

SAME  BASIC  SHAPE.   THE  SOLUTION  IS  GAINED  BY  USING  A  MATRIX  * 

OF  INFLUENCE  COEFFICIENTS  AND  THE  WING  SHAPE,  AND  SOLVING  FOR  * 

THE  VORTEX  STRENGTHS.  * 

* 

ALL  REQUIRED  PARAMETERS  AND  VARIABLES  ARE  CONTAINED  IN  COMMON.  * 
A  DESCRIPTION  OF  THE  COMMON  BLOCK  VARIABLES  IS  CONTAINED  IN  * 
OTHER  ROUTINES.  * 

* 


VARIABLES  USED  IN  THE  SUBROUTINE  ARE: 

ALF  :  ANGLE  OF  ATTACK  IN  RADIANS 

DELTA  :  THE  LEADING  EDGE  SWEEP  ANGLE  IN  RADIANS 

ADELTA  :  THE  LEADING  EDGE  SWEEP  ANGLE  IN  DEGREES 

CBAR  :  THE  AVERAGE  CHORD 

CL  :  INTERIM  VALUE,  USED  TO  COMPUTE  A  FINAL  VALUE 

CD  :  INTERIM  VALUE,  USED  TO  COMPUTE  A  FINAL  VALUE 

CM  :  INTERIM  VALUE,  USED  TO  COMPUTE  A  FINAL  VALUE 

CXJ  :  INTERIM  VALUE,  USED  TO  COMPUTE  A  FINAL  VALUE 

CZJ  :  INTERIM  VALUE,  USED  TO  COMPUTE  A  FINAL  VALUE 


* 


* 
* 

* 
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* 

* 
* 
* 
* 

* 


CMJ  :  INTERIM  VALUE,  USED  TO  COMPUTE  A  FINAL  VALUE 

XAC  :  X  POSITION  OF  THE  AERODYNAMIC  CENTER 

CLT  :  TOTAL  LIFT  COEFFICIENT  OF  THE  FLAT  WING 

CDT  :  TOTAL  DRAG  COEFFICIENT  OF  THE  FLAT  WING 

CMT  :  TOTAL  MOMENT  COEFFICIENT  OF  THE  FLAT  WING 

ABOUT  THE  AERODYNAMIC  CENTER 

CDOCL2  :  CD/CLA2,  ALSO  CALLED  K  IN  THE  DRAG  POLAR 

SUBROUTINES  CALLED  ARE: 

SET85   :  FORTRAN,  USED  TO  GENERATE  THE  WING  AND  SECTION 

ARRAY  CONTAINING  PLANFORM  COORDINATES. 
DNWS85  :  FORTRAN,  USED  TO  COMPUTE  THE  DOWNWASH  GENERATED 

AT  A  POINT  BY  A  HORSE-SHOE  VORTEX  AND  IT'S  MIRROR 

IMAGE  ON  THE  OPPOSITE  SEMI -SPAN. 
GAUSS   :  FORTRAN,  USED  TO  SOLVE  THE  INFLUENCE  COEFF  MATRIX 

AND  SLOPE  VECTOR,  FOR  THE  VORTEX  STRENGTH.   USES 

SIMPLE  GAUSS  ELIMINATION. 
MULT2   :  FORTRAN,  USED  TO  COMPUTE  THE  DOWNWASH  GENERATED 

ON  EACH  VORTEX  ELEMENT,  NECESSARY  TO  FIND  THE 

FORCES  ACTING  ON  THE  WING. 


* 
* 
* 
* 
* 
* 


* 


* 

* 
* 


SUBROUTINE  LOAD85 

REAL  LAMDA 

CHARACTER*3,  ELLIP 

COMMON  PI,B,WING(101,23) ,  SECTN(11 , 10) , AR, LAMDA, DELTA, NX, NY, 
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*  ALPHA, CLDSRD,ELLIP,XAC,CLA 
COMMON/COF/A(101 , 101) ,  NEQNS 
COMMON/CAM/AA(101,101),  AS (101 ,101) 
COMMON/FILS/IIN , IOUT , JOUT , KOUT , LOUT , MOUT , NOUT 

60     FORMAT  (  '  FLAT  WING', 

*  /'  CL     =  '.F12.6/'  CD     =  ',F12.6/'  CD/CL2  =  ',F12.6 

*  /'  CMAC    =  ',F12.6/'  XAC    =  ',F12.6/'  AR     =  ',F12.6 

*  /'  LAMDA  =  ',F12.6/'  DELTA  =  '.F12.6/'  NX     =  ',13 

*  /'NY      =  ',13    /'  ALPHA  =  '.F12.6/'  CLDSRD  =  '.F12.6 

*  /'  ELLIP  =  ',A3    /'  CLa,  Lift  Curve  Slope' 

*  /'  per/Deg=  ' ,F12.6) 
********************************************* 

*  * 

*  COMPUTE  THE  ANGLE  OF  ATTACK  IN  RADIANS.  * 

*  * 
*********************************************************************** 

ALF  =  ALPHA*PI/180.0 
SINALF  =  SIN(ALF) 
COSALF  =  COS (ALF) 
*********************************************************************** 

*  * 

*  CALL  THE  SETUP  ROUTINE.  USED  TO  BUILD  THE  COORDINATES  OF  THE 

*  WING.   THEY  ARE  STORED  IN  THE  ARRAYS  WING  AND  SECTN.  * 

*  * 
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****************************************************** 

CALL  SET 8 5 
*********************************************************************** 

*  * 

*  BUILD  THE  THREE  COEFFICIENT  MATRICES,  A(),  AA(),  AND  AS ( )  * 

*  A()  AND  AA()  ARE  THE  SAME,  AND  ARE  FOR  EVALUATING  FLOW  * 

*  TANGENCY  AT  THE  EDGE  OF  THE  GRID  ELEMENTS.  * 

*  AS()  IS  FOR  EVALUATING  DOWN-WASH  GENERATED  ON  THE  VORTEX  * 

*  ELEMENT  ITSELF. 


*********************************************************************** 
DO  130  1=1, NY 
DO  120  J=1,NX 

IJ  =  (I-1)*NX+J 
*********************************************************************** 

*  * 

*  AUGMENT  THE  A()  MATRIX  WITH  THE  RIGHT  HAND  SIDE,  WING  SHAPE.     * 

*  * 
*********************************************************************** 

A(IJ,NEQNS+1)  -  SINALF 

*  A(IJ,NEQNS+1)  =  ALF 
DO  110  K=1,NY 

DO  100  L=1,NX 

KL  =  (K-1)*NX+L 

CALL  DNWS85(IJ,KL,A(KL,IJ) , 'CONT' ) 
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AA(KL.IJ)  =  A(KL,IJ) 

CALL  DNWSSSCIJ.KL.ASCKL.IJ) , 'SELF') 

100  CONTINUE 

110  CONTINUE 

120      CONTINUE 

130    CONTINUE 

*********************************************************************** 

*  * 

*  COMPUTE  THE  AVERAGE  CHORD  * 

*********************************************************************** 

CBAR  =  (SECTN(l,l)+SECTN(NY,l))/2.0 
*********************************************************************** 

*  * 

*  CALL  THE  GAUSS  ROUTINE,  TO  SOLVE  FOR  THE  LOAD  VECTOR  * 

*  THE  RESULT  IS  RETURNED  IN  PLACE  OF  THE  ORIGINAL  LOAD  VECTOR     * 

*  * 
*********************************************************************** 

CALL  GAUSS (1) 
*********************************************************************** 

*  * 

*  STORE  THE  RESULT  IN  THE  WING()  MATRIX  * 

*  CONVERT  THE  CIRCULATION,  OR  VORTEX  STRENGTH  IN  ELEMENT  15  TO    * 

*  DELTA  CP,  IN  ELEMENT  16  * 

*  * 
*********************************************************************** 
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DO  210  I  =  1,NY 
DO  200  J  =  l.NX 

IJ  =  (I-1)*NX+J 
WING(IJ,15)  =  A(IJ,NEQNS+1) 
WING(IJ,16)  =  2.0*WING(IJ,15)/WING(IJ,13) 
200      CONTINUE 
210    CONTINUE 
************************************************** 

*  * 

*  COMPUTE  THE  DOWNWASH  ON  EACH  VORTEX  ELEMENT,  NECESSARY  TO       * 

*  FIND  THE  FORCES  ACTING  ON  THE  WING.  * 

*  * 
*********************************************************************** 

CALL  MULT 2 
*********************************************************************** 

*  * 

*  INITIALIZE  VALUES  AND  COMPUTE  THE  FORCES  ACTING  ON  THE  WING     * 

*  SECTIONS  AND  ALSO  THE  TOTAL  WING.  * 

*  * 
*********************************************************************** 

CD  =  0.0 
CL  =  0.0 
CM  =  0.0 
DO  330  I-l.NY 

CXJ  =  0.0 

CZJ  =  0.0 
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CMJ  =  0.0 

DO  320  J=1,NX 

IJ  =  (I-1)*NX+J 

CZJ  =  WING(IJ,15)*2.0*(1.-WING(IJ,22)*SINALF)/(B*CBAR) 
CXJ  =  WING(IJ,15)*WING(IJ,22)*2.0*COSALF/(B*CBAR) 
CMJ  =  -WING(IJ,15)*WING(IJ,4)*2.0* 
*  (l.-WING(IJ,22)*SINALF)/(B*CBAR*CBAR) 

SECTN(I,2)  =  SECTN(I,2)+CZJ 
SECTN(I,3)  =  SECTN(I,3)+CXJ 
SECTN(I,4)  =  SECTN(I,4)+CMJ 
320      CONTINUE 

CL  =  CL+SECTN(I,2)*SECTN(I,6) 
CD  =  CD+SECTN(I,3)*SECTN(I,6) 
CM  =  CM+SECTN(I,4)*SECTN(I,6) 
SECTN(I,5)  -  -SECTN(I,4)/SECTN(I,2) 
330    CONTINUE 
*********************************************************************** 

*  * 

*  NOW  COMPUTE  THE  MOMENTS  ABOUT  THE  AERODYNAMIC  CENTER  * 

*  * 
*********************************************************** 

XAC  =  -CM*CBAR/CL 

DO  335  I  =  1,NY 

SECTN(I,4)  =  SECTN(I,4)  +  SECTN(I , 2)*XAC/CBAR 
335    CONTINUE 
*********************************************************************** 
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*  COMPUTE  THE  FINAL  TOTAL  WING  COEFFICIENTS  AND  WRITE  TO  THE      * 

*  FILE.  * 

*  * 

CLT  =  2.0*CL 
CDT  =  2.0*CD 

CMT  =  2.0*CM  +  (CLT*XAC/CBAR) 
CLA  =  CLT/ALPHA 
CDOCL2  =  CDT/CLT**2 
ADELTA  =  DELTA  *  180. /PI 
OPEN  (UNIT  =  MOUT,  FILE  =' FLAT. DAT') 

WRITE (MOUT , 60 )  CLT , CDT , CDOCL2 , CMT , XAC , AR , LAMDA , ADELTA , NX , NY , 
*  ALPHA,  CLDSRD,  ELLIP,  CLA 
CLOSE  (MOUT) 
340   RETURN 
END 
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SHAP85 


SUBROUTINE  SHAP85 


SUBROUTINE  TO  DEVELOP  A  WING  SHAPE  THAT  IS  ASSOCIATED  WITH  * 

A  SPECIFIED  LOADING  FUNCTION.   THE  LOADING  FUNCTION  IS  USUALLY  * 

ELLIPTIC,  IN  THE  SPAN  AND  CHORD  DIRECTION.  THE  SHAPE  WILL  * 

USUALLY  RESEMBLE  A  PARABOLIC  CAMBER,  IN  CHORD,  AND  HAVE  SOME  * 

DEGREE  OF  TWIST,  ALONG  THE  SPAN.  * 

* 

ALL  REQUIRED  PARAMETERS  AND  VARIABLES  ARE  CONTAINED  IN  COMMON.  * 

A  DESCRIPTION  OF  THE  COMMON  BLOCK  VARIABLES  IS  CONTAINED  IN  * 

OTHER  ROUTINES.  * 

* 

VARIABLES  USED  IN  THE  SUBROUTINE  ARE:  * 

CBAR    :  THE  AVERAGE  CHORD  * 

CLILST  :  CLREF  FROM  THE  PREVIOUS  ITERATION 

CLREF     :  A  SCALER,  USED  TO  INCREASE  OR  DECREASE  THE  AMPL  * 

OF  THE  SHAPE,  IN  ORDER  TO  ACHIEVE  THE  DESIRED  * 

LIFT  COEFFICIENT.  * 

MOD     :  A  COUNTER  USED  TO  LIMIT  THE  NUMBER  OF  ITERATIONS  * 

THAT  THE  PROGRAM  CAN  USE  IN  AN  ATTEMPT  TO  GENERATE  * 

A  DESIRED  LIFT  COEFFICIENT.  * 

H       :  A  TEMPORARY  VARIABLE  USED  WHILE  INTEGRATING  THE 

HEIGHT  OF  THE  WING  ALONG  A  CHORD.  * 
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* 

* 


COSCAM 


SINCAM 


THE  COSINE  OF  THE  ANGLE  MADE  BY  AN  INDIVIDUAL  GRID   * 

ELEMENT  AND  THE  REMOTE  VELOCITY.  * 

THE  SIN  OF  THE  ANGLE  MADE  BY  AN  INDIVIDUAL  GRID  * 

ELEMENT  AND  THE  REMOTE  VELOCITY.  * 

INTERIM  VALUE,  USED  TO  COMPUTE  A  FINAL  VALUE  * 

INTERIM  VALUE,  USED  TO  COMPUTE  A  FINAL  VALUE  * 

INTERIM  VALUE,  USED  TO  COMPUTE  A  FINAL  VALUE  * 

X  POSITION  OF  THE  CENTER  OF  PRESSURE  * 

TOTAL  LIFT  COEFFICIENT  OF  THE  CAMBERED  WING.  * 

TOTAL  DRAG  COEFFICIENT  OF  THE  CAMBERED  WING.  * 

TOTAL  MOMENT  COEFFICIENT  OF  THE  CAMBERED  WING.  * 

ABOUT  THE  AERODYNAMIC  CENTER  * 

CD/CLA2,  ALSO  CALLED  K  IN  THE  DRAG  POLAR  * 

SUBROUTINES  CALLED  ARE:  * 

MULT1   :  FORTRAN,  USED  TO  COMPUTE  THE  SHAPE  OF  THE  WING  * 

FROM  THE  LOADING  VECTOR  AND  THE  INFLUENCE  COEFF  * 

MATRIX  * 

MULT3   :  FORTRAN,  USED  TO  COMPUTE  THE  DOWNWASH  GENERATED  * 

ON  EACH  VORTEX  ELEMENT,  NECESSARY  TO  FIND  THE  * 

FORCES  ACTING  ON  THE  WING.  * 


* 

CCL 

* 

CCD 

* 

CCM 

* 

CXCP 

* 

CCLT 

* 

CCDT 

* 

CCMT 

* 

* 

CCDOCL 

SUBROUTINE  SHAP85 
REAL  LAMDA 
CHARACTER*3,  ELLIP 


123 


COMMON  PI,B,WING(101,23)  ,  SECTN(11 , 10) , AR, LAMDA, DELTA, NX, NY, 
*  ALPHA ,CLDSRD,ELLIP,XAC,CLA 
COMMON/COF/A(101,101) ,  NEQNS 
COMMON/CAM/ AA( 10 1,101) ,  AS (101, 101) 
COMMON/FILS/I IN , IOUT , JOUT , KOUT , LOUT , MOUT , NOUT 
CLILST  =0.0 
CLREF  =0.05 
MOD  =  0 


66 

FORMAT  ('  UN 

ABLE 

70 

FORMAT  (  '  C 

AMBE 

* 

/' 

CL 

=  ' 

* 

/' 

CMAC 

=  ' 

* 

/' 

LAMDA 

=  ' 

* 

/' 

NY 

=  ' 

* 

/' 

ELLIP 

=  ' 

,F12.6/'  CD     =  '.F12.6/'  CD/CL2  =  '.F12.6 

.F12.6/  '  AR     =  ' ,F12.6 

.F12.6/'  DELTA  =  '.F12.6/'  NX     =  ',13 

,13    /'  ALPHA  =  '.F12.6/'  CLDSRD  =  '.F12.6 
,A3) 

*********************************************************************** 

*  * 

*  COMPUTE  THE  AVERAGE  CHORD  * 

*  * 
*********************************************************************** 

CBAR  =  (SECTN(l,l)+SECTN(NY,l))/2.0 
******************************************************** 

*  * 

*  USE  THE  SCALER,  CLREF,  TO  PUT  AN  ELLIPTIC  LOAD  DISTRIBUTION 

*  IN  THE  VECTOR  (I J, 17).   AT  THE  SAME  TIME  COMPUTE  THE  DELTA  CP    * 

*  IN  THE  VECTOR  (I J, 18).  * 


1" 


****************************************************** 
390    DO  400  I J  =  l.NEQNS 

WING(IJ,17)  =  WING(IJ,21)*CLREF 

WING(IJ,18)  =  2.0*WING(IJ,17)/WING(IJ,13) 
400    CONTINUE 
*********************************************************************** 

*  * 

*  CALL  THE  SUBROUTINE  MULT1  WHICH  WILL  * 

*  MULTIPLY  THE  LOAD  VECTOR,  (I J, 17)  BY  THE  COEFFICIENT  MATRIX,  * 

*  AA(  ),  TO  GET  THE  SHAPE  VECTOR,  (IJ.19)  * 

*  * 
*********************************************************************** 

CALL  MULT1 
*********************************************************************** 

*  * 

*  CHECK  THE  RESULTING  SHAPE  VECTOR  (I J, 19),  TO  SEE  IF  THE  ANGLE  * 

*  THAT  WAS  REQUIRED  HAS  A  SIN  GREATER  THAN  1,  THAT  IS  SIMPLY  * 

*  CHECKING  THE  SHAPE  VECTOR,  TO  SEE  IF  IT  IS  GREATER  THAN  1.  * 

*  IF  IT  IS,  THE  SCALER  MULTIPLIER,  CLREF,  IS  DECREASED.   AT  THE 

*  SAME  TIME,  MOD,  THE  ITERATION  LIMITER,  IS  INCREMENTED.   IF  * 

*  MOD  IS  GREATER  THAN  5,  THAT  MEANS  WE  OVERSHOT  5  TIMES,  AND  * 

*  PROBABLY  WILL  NEVER  GET  TO  THE  DESIRED  LIFT  COEFFICIENT.  * 

*  IN  THAT  CASE,  COMPUTE  THE  LAST  VALID  SHAPE,  AND  PROVIDE  THOSE  * 

*  DATA .  * 

*  * 
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*********************************************************************** 

DO  410  IJ  =  1, NEONS 

IF  (ABS(WING(IJ,19)) .GT.1.0)  THEN 
CLREF  ■=  CLREF*.9 
MOD  =  MOD  +  1 
IF  (MOD.GT.5)  THEN 
WRITE  (IOUT.66) 
CALL  DELAY(6) 
GO  TO  540 
ELSE 

GO  TO  390 
END  IF 
ELSE 
ENDIF 
410    CONTINUE 
****************************************** 

*  * 

*  BEGIN  TO  COMPUTE  THE  FORCE  AND  MOMENT  COEFFICIENTS  ON  THE  * 

*  CAMBERED  WING.   MULT3  WILL  MULTIPLY  THE  SPECIFIED  LOADING  * 

*  BY  THE  COEFFICIENT  MATRIX,  AS,  TO  GET  THE  DOWNWASH  AT  THE  * 

*  CENTER  OF  EACH  CAMBERED  ELEMENT.   INTERMEDIATE  VARIABLES  ARE 

*  INITIALIZED.  * 

*  COEFFICIENTS  FOR  SECTIONS  AND  THE  WHOLE  WING  ARE  COMPUTED.  * 

*  * 
*********************************************************************** 

CALL  MULT 3 
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CCD  =  0.0 
CCL  =0.0 
CCM  =0.0 
DO  530  1=1, NY 
CCXJ  =0.0 
CCZJ  =0.0 
CCMJ  =0.0 
SECTN(I,7)  =  0.0 
SECTN(I,8)  =  0.0 
SECTN(I,9)  =  0.0 
H  =  0.0 

IJ  =  (I-1)*NX+1 
DO  520  J=1,NX 

IJ  =  (I-1)*NX+J 
H  =  H-WING(IJ,19)*WING(IJ,13) 
WING(IJ,20)  =  H 

COSCAM  =  COS(ASIN(WING(IJ,19))) 
SINCAM  =  WING(IJ,19) 

CCZJ  =  WING(IJ,17)*2.0*(1.-WING(IJ,23)*SINCAM)/(B*CBAR) 
CCXJ  =  WING(IJ,17)*WING(IJ,23)*2.0*COSCAM/(B*CBAR) 
CCMJ  =  -WING(IJ,17)*WING(IJ,4)*2.0* 
*  (1. -WING(IJ,23)*SINCAM)/(B*CBAR*CBAR) 

SECTN(I,7)  =  SECTN(I,7)+CCZJ 
SECTN(I,8)  =  SECTN(I,8)+CCXJ 
SECTN(I,9)  =  SECTN(I,9)+CCMJ 
520      CONTINUE 
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CCL  =  CCL+SECTN(I,7)*SECTN(I,6) 
CCD  =  CCD+SECTN(I,8)*SECTN(I,6) 
CCM  -  CCM+SECTN(I,9)*SECTN(I,6) 
*********************************************************************** 

*  * 

*  COMPUTE  THE  SECTION  TWIST.   THE  CAMBERED  WING  CL  * 

*  WILL  BE  THE  SAME  AS  CLDSRD ,  WITHIN  A  VERY  SMALL  MARGIN.         * 

*  * 

*********************************************************************** 

SECTN(I.IO)  =  -WING(I*NX,20)*180./(SECTN(I,1)*PI) 
********************************************** 

*  * 

*  COMPUTE  THE  MOMENT  COEFFICIENT  ABOUT  THE  AERODYNAMIC  CENTER.     * 

*  * 
*********************************************************************** 

SECTN(I,9)  =  SECTN(I,9)  +  (SECTN(I , 7)*XAC/CBAR) 
530    CONTINUE 

CXCP  =  -CCM/CCL 

CCLT  =  2.0*CCL 

CCDT  =  2.0*CCD 

CCMT  =  2.0*CCM  +  (CCLT*XAC/CBAR) 

CCDOCL  =  CCDT/CCLT**2 

CLILST  =  CLREF 

IF  (ABS (CLDSRD- CCLT) . LT . (CLDSRD*0 . 01) )  GO  TO  540 

CLREF  =  CLREF*CLDSRD/CCLT 

GO  TO  390 
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******************************************************* 

*  * 

*  THE  DESIRED  LIFT  COEFFICIENT  HAS  BEEN  ACHIEVED,  WITHIN  A  * 

*  A  TOLERANCE  OF  +/"  1.0%.  THE  CORRECT  VORTEX  STRENGTH  AND  * 

*  DELTA  CP  ARE  PUT  INTO  THE  WING  ARRAY,  AND  THE  FINAL  CAMBERED  * 

*  WING  VALUED  OF  LIFT,  DRAG,  K,  MOMENT,  AND  CENTER  OF  PRESSURE  * 

*  ARE  PRINTED.  * 

*  * 
*********************************************************************** 
540    CONTINUE 

DO  550  I J  =  l.NEQNS 

WING(IJ,17)  =  WING(IJ,21)*CLILST 

WING (IJ, 18)  =  2.0*WING(IJ,17)/WING(IJ,13) 
550    CONTINUE 
*********************************************************************** 

*  * 

*  WRITE  THE  TOTAL  WING  COEFFICIENTS  TO  THE  SCREEN  * 

*  * 

***************************************************** ** * *************** 
OPEN  (UNIT  =  MOUT,  FILE  =  'CAMBER.DAT') 
ADELTA  =  DELTA  *180./PI 

WRITE (MOUT, 70)  CCLT , CCDT, CCDOCL, CCMT ,      AR , LAMDA , ADELTA , NX , NY , 
*  ALPHA,  CLDSRD,  ELLIP 
CLOSE  (MOUT) 
RETURN 
END 
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SET85 

************************************************* 

*  * 

*  SUBROUTINE  SET85  * 

*  * 

*  This  program  is  written  for  the  setup  of  a  straighthorse  shoe    * 

*  vortex  over   tapered  and  swept  wing.  It  uses  cosine  spacing  * 

*  for  both  the  span  and  chord  elements.   This  concentrates  * 

*  elements  near  the  edges  of  the  wing,  where  they  are  most  * 

*  important.  * 

*  * 

*  Column  elements  in  the  WING  array  are:  * 

*  1   sequential  position,  also  called  IJ  * 

*  2   integer  y  position  * 

*  3   integer  x  position  * 

*  4   x  center  of  element  * 

*  5   y  center  of  element  * 

*  6   x  position  of  horse  shoe  vortex  * 

*  7   blank  * 

*  8   ya  position  of  one  corner  of  horse  shoe  vortex  * 

*  9   yb  position  of  one  corner  of  horse  shoe  vortex  * 

*  10  xp  position  of  flow  tangency  point,  3/4  chord  * 

*  11  xp  position  of  down  wash  point,  1/4  chord  * 

*  12  yp  position  of  flow  tangency  or  down  wash  pt,  center  of  elem  * 

*  13  del  x,  chord  wise  dimension  * 

*  14  del  y,  span  wise  dimension  * 
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*  15  flat  plate  vorticity  * 

*  16  flat  plate  delta  cp  * 

*  17  cambered  voricity,  adjusted  for  CLi  * 

*  18  cambered  delta  cp  * 

*  19  cambered  slope  * 

*  20  cambered  height  * 

*  21  cambered  vorticity,  if  CLREF  were  =1.0  * 

*  22  downwash  at  center  of  flat  plate  element  * 

*  23  downwash  at  center  of  cambered  element  * 

*  * 

*  Column  elements  in  the  SECTN  array  are:  * 

*  1   chord  * 

*  2   flat  wing  lift  force  per  unit  span  * 

*  3   flat  wing  drag  force  per  unit  span  * 

*  4   flat  wing  moment  about  the  aerodynamic  center  * 

*  5   flat  wing  XCP,  center  of  pressure  relative  to  x=0  * 

*  6   section  delta  y  * 

*  7   cambered  wing  lift  force  per  unit  span  * 

8  cambered  wing  drag  force  per  unit  span  * 

9  cambered  wing  moment  about  the  aerodynamic  center  * 

10  cambered  wing  twist.  * 


VARIABLES  USED  ARE:  * 

B       :   WING  SPAN  * 

DTHE    :   AN  ANGULAR  ELEMENT  WIDTH,  IN  RADIANS   SPAN          * 

DPSI    :   AN  ANGULAR  ELEMENT  DEPTH,  IN  RADIANS,  CHORD         * 


* 
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* 
* 

* 
* 
* 

* 
* 
* 

* 
* 


THEA  :  INTERMEDIATE  VARIABLE 
THEA  :  INTERMEDIATE  VARIABLE 
PS I A  :  INTERMEDIATE  VARIABLE 
PS IB  :  INTERMEDIATE  VARIABLE 
ELE  :  INTERMEDIATE  VARIABLE 
ELT  :  INTERMEDIATE  VARIABLE 
PSI     :   CHORDWISE  ANGULAR  COORD  FOR  THE  CENTER  OF  AN 

ELEMENT.   USED  IN  GENERATING  ELLIPTIC  LOAD. 
THETA   :   SPANWISE  ANGULAR  COORD  FOR  THE  CENTER  OF  AN 

ELEMENT.  USED  IN  GENERATING  ELLIPTIC  LOAD. 

FUNCTIONS  USED  ARE: 

BETA    :   COMPUTES  THE  LENGTH  OF  A  LOCAL  SECTION  CHORD 

ALL  VARIABLES  AND  PARAMETERS  ARE  PASSED  IN  COMMON  BLOCKS 
THE  ONLY  VARIABLES  CHANGED  IN  THIS  ROUTINE  ARE: 

B 

WING() 

SECTNQ 


* 

* 
* 
* 

* 
* 


* 
* 
* 
* 


************************************************ 

SUBROUTINE  SET85 

COMMON  PI,B,WING(101,23) ,  SECTN(11 , 10) , AR.LAMDA.DELTA.NX.NY, 
*  ALPHA , CLDSRD , ELLI P , XAC , CLA 

REAL  LAMDA 
*********************************************************************** 
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*  * 

*  DEFINE  THE  FUNCTION  BETA,  WHICH  COMPUTES  THE  LENGTH  OF  A  LOCAL  * 

*  CHORD .  * 

************************************************ 

BETA(Y,LAMDA,AR)  =  1 . 0- (Y*4 . 0*(1 . 0-LAMDA)/(AR*(l . 0+LAMDA) ) ) 
*********************************************************************** 

*  * 

*  DEFINE  THE  SIZE  OF  SPAN  AND  CHORD  WISE  GRID  ELEMENTS.  * 

*  DETERMINE  THE  TOTAL  WING  SPAN.  * 

*  * 
*********************************************************************** 

*  DTHE  =  PI*0.5/FLOAT(NY) 
*********************************************************************** 

*  * 

*  THIS  STATEMENT  SETS  THE  GRID  UP  SO  THAT  THE  OUTER  ELEMENT       * 

*  IS  IGNORED.  * 

*  * 
****************************************************** ***************** 

DTHE  =  PI*0.5/FLOAT(NY+1) 
DPS I  =  PI/FLOAT (NX) 
B  =  AR*(1.0+LAMDA)/2.0 
DO  110  I  =  1,NY 
DO  100  J  =  1,NX 

*********************************************************************** 
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*  COMPUTE  ALL  REQUIRED  COORDINATES  FOR  THE  WING()  ARRAY.  * 

*  * 
************************************************ 

IJ  =  (I-1)*NX+J 

WING(IJ.l)  =  FLOAT(IJ) 

WING(IJ,2)  =  FLOAT(I) 

WING (I J, 3)  =  FLOAT(J) 

THEA  =  (FLOAT(I-1)*DTHE)+(PI/2.0) 

THEB  =  (FLOAT(I)*DTHE)+(PI/2.0) 

PSIA  =  (FLOAT(J-l)*DPSI) 

PSIB  =  (FLOAT(J)*DPSI) 

WING(IJ,8)  =  -B*COS(THEA)/2.0 

WING (I J, 9)  =  -B*COS(THEB)/2.0 

WING(IJ,12)  =  -B*(COS(THEA)+COS(THEB))/4.0 

ELE  =  (1.0  -  COS(PSIA))* 

*  BETA (WING (I J, 12) , LAMDA, AR)/2 . 0 
ELT  =  (1.0  -  COS(PSIB))* 

*  BETA(WING(IJ,12) , LAMDA, AR)/2.0 

*  THESE  TWO  STATEMENTS  SET  THE  POINTS  AT  1/4  AND  3/4  CHORD 

WING(IJ.IO)  =  WING(IJ, 12) *TAN (DELTA) +(ELE+3*ELT)/4 . 0 

WING(IJ.ll)  =  WING(IJ,12)*TAN(DELTA)+(3*ELE+ELT)/4.0 

WING(IJ,4)  =  WING(IJ,12)*TAN(DELTA)+(ELE+ELT)/2.0 

WING(IJ,6)  =  WING(IJ,11) 

WING(IJ,5)  =  WING(IJ,12) 

WING(IJ,13)  =  ELT-ELE 

WING(IJ,14)  =  WING(IJ,9)-WING(IJ,8) 

134 


PSI  =  ACOS(1.0-(2.0*(WING(IJ,H)-(WING(IJ,5)* 

*  TAN ( DELTA))) /BETA (WING (IJ, 5) ,LAMDA,AR))) 
THETA  =  ACOS(-2.0*WING(IJ,5)/B) 

WING(IJ,21)  =  4.0*(1.0+LAMDA)*SIN(THETA)*SIN(PSI)* 

*  WING(IJ,13)/(PI**2*BETA(WING(IJ,5) ,LAMDA,AR)) 
WING(IJ,7)  =0.0 

WING(IJ,15)  =  0.0 

WING (I J, 16)  =0.0 

WING(IJ,17)  =  0.0 

WING (IJ, 18)  =0.0 

WING (I J, 19)  =0.0 

WING(IJ,20)  =  0.0 

WING (IJ, 22)  =0.0 

WING(IJ,23)  =  0.0 
100      CONTINUE 
************************************************ 

*  * 

*  COMPUTE  ALL  REQUIRED  COORDINATES  FOR  THE  SECTN()  ARRAY.         * 

*  * 
***********************************************  ******************************* 

SECTN(I.l)  =  BETA(WING(IJ,12) ,LAMDA,AR) 

SECTN(I,2)  =  0.0 

SECTN(I,3)  =  0.0 

SECTN(I,4)  =  0.0 

SECTN(I,5)  =0.0 

SECTN(I,6)  =  WING(IJ,14) 
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SECTN(I,7)  =  0.0 
SECTN(I,8)  =  0.0 
SECTN(I,9)  =  0.0 
SECTN(I.IO)  =  0.0 
110    CONTINUE 

RETURN 

END 


1" 


DNWS85 

*************************************************** 

*  * 

*  SUBROUTINE   DNWS85 (IJ ,KL,W, IND)  * 

*  * 

*  COMPUTE  DOWNWASH  ON  GRID  ELEMENT  NUMBER  KL  DUE  TO  * 

*  VORTICES  ON  GRID  ELEMENT  IJ,  AND  IT'S  MIRROR  IMAGE  * 

*  * 

*  VARIABLES  USED  ARE:  * 

*  I J      :   VORTEX  GRID  ELEMENT  NUMBER,  INPUT  * 

*  KL      :   DOWNWASH  GRID  ELEMENT  NUMBER,  INPUT  * 

*  W       :   THE  DOWN  WASH,  RETURNED  TO  THE  CALLING  PROGRAM  * 

*  OUTPUT  * 

*  IND     :   A  FLAG  TO  DETERMINE  WHETHER  TO  USE  TRAILING  * 

*  EDGE  OF  THE  GRID  ELEMENT,  OR  THE  CENTER  OF  * 

*  THE  ELEMENT,  INPUT  * 

*  * 

*  FUNCTIONS  USED  ARE:  * 

*  WHV1    :   FORTRAN,  COMPUTES  THE  DOWNWASH  DUE  TO  ONE  CORNER   * 

*  OF  THE  HORSESHOE  VORTEX  * 

*  WHV2    :   FORTRAN,  COMPUTES  THE  DOWNWASH  DUE  TO  JUST  THE  * 

*  TRAILING  FILAMENTS  OF  THE  HORSESHOE  VORTEX  * 

*  * 

*********************************************************************** 

SUBROUTINE  DNWS85(IJ ,KL,W, IND) 
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COMMON  PI,B,WING(101,23),  SECTN(11 , 10) , AR.LAMDA, DELTA.NX.NY, 

*  ALPHA, CLDSRD,ELLIP,XAC,CLA 
CHARACTER*4  IND 
IF(IND.EQ.'CONT')  GOTO  100 
IF(IND.EQ. 'SELF')  GOTO  110 

100   W  =  WHV1(WING(KL,10) ,WING(KL,12) ,WING(IJ,6) ,WING(IJ,8)) 

*  -  WHV1(WING(KL,10) ,WING(KL,12) ,WING(IJ,6) ,WING(IJ,9)) 

*  -  WHV1(WING(KL,10) ,WING(KL,12) ,WING(IJ,6) , -WING(IJ,8)) 

*  +  WHV1(WING(KL,10) ,WING(KL,12) ,WING(IJ,6) , -WING(IJ,9)) 
W  =  W*.25/PI 

RETURN 
110   W  =  WHV2(WING(KL,11) ,WING(KL,12) ,WING(IJ , 6) ,WING(IJ , 8) ) 

*  -  WHV2(WING(KL,11) ,WING(KL,12) ,WING(IJ,6) ,WING(IJ,9)) 

*  -  WHV2(WING(KL,11) ,WING(KL,12) ,WING(IJ,6) , -WING(IJ , 8) ) 

*  +  WHV2(WING(KL,11),WING(KL,12),WING(IJ,6),-WING(IJ,9)) 
120   W  =  W*.25/PI 

RETURN 
END 
******************************************************* 

*  * 

*  FUNCTION  WHV1(X1,Y1,X2,Y2)  * 

*  * 
*********************************************************************** 

FUNCTION  WHV1(X1,Y1,X2,Y2) 

IF  (ABS(X1-X2).LT. .0001)  GOTO  100 

WHV1  =  (1.0+SQRT((X1-X2)**2  +  (Y1-Y2)**2)/(X1-X2))/(Y1-Y2) 
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RETURN 
100   WHV1  -  1.0/(Y1-Y2) 
RETURN 
END 

*********************************************************************** 

*  FUNCTION  WHV2(X1,Y1,X2,Y2)  * 

*  * 
************************************************* 

FUNCTION  WHV2(X1,Y1,X2,Y2) 
IF  (ABS(X1-X2).LT. .0001)  GOTO  100 

WHV2  =  (1.0+(X1-X2)/SQRT((X1-X2)**2  +  (Yl-Y2)**2) )/(Yl-Y2) 
RETURN 
100   WHV2  =  1.0/(Y1-Y2) 
RETURN 
END 
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MULT1 

*********************************************** 

*  * 

*  SUBROUTINE  MULT1  * 

*  * 

*  SUBROUTINE  TO  MULTIPLY  THE  SPECIFIED  LOADING  VECTOR,  (1,17)      * 

*  BY  THE  COEFFICIENT  MATRIX,  AA(  ),  TO  GET  THE  WING  SHAPE         * 

*  VECTOR  (1,19)  .  * 

*  * 

*  VARIABLE  NAMES  ARE  CONTAINED  IN  THE  MAIN  PROGRAM.  * 

*  ALL  VARIABLES  AND  PARAMETERS  ARE  PASSED  IN  COMMON.  * 

*  THE  ONLY  VARIABLE  CHANGED  IS  WING (I, 19)  * 

*  * 
*********************************************************************** 

SUBROUTINE  MULT1 
REAL  LAMDA 
CHARACTER* 3  ELLIP 

COMMON  PI,B,WING(101,23) ,  SECTN(11 , 10) , AR, LAMDA, DELTA, NX, NY, 
*  ALPHA, CLDSRD,ELLIP,XAC,CLA 
COMMON/COF/A(101 , 101) ,  NEQNS 
COMMON/CAM/AA(101,101) ,  AS(101,101) 
DO  110  I  =  1, NEQNS 
WING (I, 19)  =0.0 
DO  100  J  =  1, NEQNS 

WING(I,19)  =  WING(I,19)+AA(I,J)*WING(J,17) 
100      CONTINUE 
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110    CONTINUE 
RETURN 
END 
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MULT2 

********************************************************************** 

*  * 

*  SUBROUTINE  MULT2  * 

*  * 

*  SUBROUTINE  TO  MULTIPLY  THE  FLAT  PLATE  LOAD  VECTOR,  (1,15)      * 

*  BY  THE  COEFFICIENT  MATRIX,  AS (  ),  TO  GET  THE  DOWN  WASH         * 

*  VECTOR  (1,22)  .  * 

*  * 

*  VARIABLE  NAMES  ARE  CONTAINED  IN  THE  MAIN  PROGRAM.  * 

*  ALL  VARIABLES  AND  PARAMETERS  ARE  PASSED  IN  COMMON.  * 

*  THE  ONLY  VARIABLE  CHANGED  IS  WING (I, 22)  * 

*  * 
********************************************* 

SUBROUTINE  MULT2 
REAL  LAMDA 
CHARACTER*3,ELLIP 

COMMON  PI,B,WING(101,23) ,  SECTN(11 , 10) , AR, LAMDA, DELTA, NX, NY, 
*  ALPHA, CLDSRD,ELLIP,XAC, CIA 
COMMON/COF/A(101 , 101) ,  NEQNS 
COMMON/CAM/AA(101,101) ,  AS (101, 101) 
DO  110  I  =  1, NEQNS 
WING(I,22)  =  0.0 
DO  100  J  =  1, NEQNS 

WING(I,22)  -  WING(I,22)+AS(I,J)*WING(J,15) 
100       CONTINUE 


y- 


110    CONTINUE 
RETURN 
END 
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MULT3 

*********************************************************************** 

*  * 

*  SUBROUTINE  MULT 3  * 

*  * 

*  SUBROUTINE  TO  MULTIPLY  THE  SPECIFIED  LOADING  VECTOR,  (1,17)      * 

*  BY  THE  COEFFICIENT  MATRIX,  AS (  ),  TO  GET  THE  DOWN  WASH  * 

*  VECTOR  (1,23) .  * 

*  * 

*  VARIABLE  NAMES  ARE  CONTAINED  IN  THE  MAIN  PROGRAM.  * 

*  ALL  VARIABLES  AND  PARAMETERS  ARE  PASSED  IN  COMMON.  * 

*  THE  ONLY  VARIABLE  CHANGED  IS  WING (I, 23)  * 

*  * 
************************************************ 

SUBROUTINE  MULT3 
REAL  LAMDA 
CHARACTER* 3  ELLIP 

COMMON  PI,B,WING(101,23) ,  SECTN(11 , 10) , AR, LAMDA, DELTA, NX, NY, 
*  ALPHA, CLDSRD,ELLIP,XAC,CLA 
COMMON/COF/A( 10 1,101) ,  NEQNS 
COMMON/CAM/AA(101,101) ,  AS(101,101) 
DO  110  I  =  1, NEQNS 
WING (I, 23)  =0.0 
DO  100  J  =  1, NEQNS 

WING(I,23)  -  WING(I,23)+AS(I,J)*WING(J,17) 
100      CONTINUE 
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110   CONTINUE 
RETURN 
END 
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GRAP85 


* 
* 
* 
* 
* 
■k 
* 


PROGRAM  GRAP85 


* 
* 


THIS  IS  A  PROGRAM  WHICH  GENERATES  A  MENU  OF  GRAPHIC  OPTIONS  * 
AND  WRITES  FILES  TO  DISK  WHICH  ARE  USED  BY  GRAPHER,  A  SOFTWARE  * 
PROGRAM  FROM  GOLDEN  SOFTWARE,  INC.  BOX  281,  GOLDEN,  COLO.,  * 
80402.  IF  YOU  DON'T  HAVE  A  COPY  OF  THAT  SOFTWARE,  YOU  MUST  * 
USE  SOME  OTHER  METHOD  OF  PRESENTING  THE  DATA.  * 

* 

DATA  FILES  USED  FOR  GRAPHIC  OUTPUT  ARE:  * 

30  =  X,Y  DATA  FILE  * 

32  =  GRAPH  FILE  * 

* 

* 

ALL  PARAMETERS  ARE  PASSED  IN  COMMON  BLOCKS.  * 

VARIABLES  USED  ARE:  * 


ELLIP  :  'YES'  OR  'NO  ',  DEPENDING  ON  WHETHER  THE  SHAPE     * 

WAS  COMPUTED  * 

TYPE  :  CHARACTER  STRING  'FLAT     '  OR  'CAMBERED'          * 

ITYPE  :  INTEGER,  1  IF  FLAT,  2  IF  CAMBERED                  * 

TITLE  :  THE  TITLE  OF  THE  GRAPH,  6  POSSIBLE  TITLES  ARE  USED  * 

ISELEC  :  INTEGER  RESPONSE  FROM  THE  GRAPHIC  SCREEN           * 

XMIN  :  MIN  VALUE  ON  HORIZONTAL  AXIS  FOR  PLOT  OF  A  SECTION  * 
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* 

XMAX 

* 

YMIN 

* 

YMAX 

* 

VMAX 

* 

VMIN 

* 

VERT 

* 

IOFFS 

* 

ISECT 

* 

ESC 

* 

HSTART 

* 

* 

VINC 

* 

XINC 

* 

YINC 

* 

PLE 

* 

PTE 

* 

XLE 

* 

XTE 

* 

AXIS 

TITLE 


TITLE2 


TYPE 


MAX  VALUE  ON  HORIZONTAL  AXIS  FOR  PLOT  OF  A  SECTION  * 

MIN  VALUE  ON  HORIZONTAL  AXIS  FOR  PLOT  OF  A  SPAN  * 

MAX  VALUE  ON  HORIZONTAL  AXIS  FOR  PLOT  OF  A  SPAN  * 

MAX  VALUE  ON  VERTICAL  AXIS  OF  PLOT  * 

MIN  VALUE  ON  VERTICAL  AXIS  OF  PLOT  * 

A  TEMPORARY  VARIABLE  USED  TO  HOLD  A  VALUE.  * 
AN  OFFSET  USED  TO  LOCATE  THE  CORRECT  COLUMN  IN  THE  * 

WING  MATRIX  * 

SPANWISE  SECTION  TO  BE  PLOTTED  * 

THE  ESCAPE  CHARACTER  (27)  * 

POSITION  ON  THE  GRAPH  WHERE  THE  HORIZONAL  AXIS  * 

STARTS  * 

VERTICLE  AXIS  INCREMENT  IN  UNITS  * 

HORIZONTAL  AXIS  INCREMENT  WHEN  PLOTTING  CHORD  * 

HORIZONTAL  AXIS  INCREMENT  WHEN  PLOTTING  SPAN  * 

PLOT  LEADING  EDGE  POSITION  IN  INCHES  * 

PLOT  TRAILING  EDGE  POSITION  IN  INCHES  * 

LEADING  EDGE  POSITION  IN  PLANFORM  UNITS  * 

TRAILING  EDGE  POSITION  IN  PLANFORM  UNITS  * 

CHARACTER  ARRAY  WHICH  HOLDS  THE  TITLES  FOR  THE  * 

HORIZONTAL  AXIS.  * 

CHARACTER  ARRAY  WHICH  HOLDS  THE  TITLES  FOR  THE  * 

VERTICAL  AXIS.  * 

CHARACTER  ARRAY  WHICH  HOLDS  PART  OF  THE  GRAPH  * 

TITLE.  * 

CHARACTER  ARRAY  WHICH  HOLDS  PART  OF  THE  GRAPH  * 
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* 

* 
* 

* 
* 
* 
* 


TITLE  * 

COEF    :   CHARACTER  ARRAY  WHICH  HOLDS  PART  OF  THE  GRAPH  * 

TITLE.  * 

COEFF   :   REAL  ARRAY  WHICH  HOLDS  THE  LIFT,  DRAG  OR  MOMENT  * 

ASSOCIATED  WITH  A  PARTICULAR  GRAPH.  * 

* 

SUBROUTINES  CALLED  ARE:  * 

CDAT85  ,   FORTRAN,  GENERAL  DATA  ENTRY  * 

DELAY   ,   FORTRAN,  SLOW  DOWN  THE  PROGRAM  SO  ERROR  MESSAGES  * 

CAN  BE  SEEN  BEFORE  BEING  ERASED.  ' 


*  SDIG    ,   FORTRAN,  USED  TO  FIND  THE  NUMBER  OF  SIGNIFICANT    * 

*  DIGITS  IN  THE  AXIS  LIMITS  AND  MAKE  AXIS  LIMITS  * 

*  MORE  UNIFORM.  * 

*  RNDUP   ,   FORTRAN,  USED  TO  ROUND  THE  AXIS  LIMITS  TO  THE  NEXT  * 

*  WHOLE  DIGIT.  * 

*  * 

PROGRAM  GRAP85 
REAL  LAMDA 
CHARACTER*3  ELLIP 
CHARACTER*8  TYPE(2) 
CHARACTER*20  TITLE(6) 
CHARACTER*8  TITLE2(2,2) 
CHARACTER*!  ESC 


1' 


CHARACTER*23  AXIS(2) 
CHARACTER* 7  COEF(2,3) 
REAL  C0EFF(2,4) 

COMMON  PI,B,WING(101,23),  SECTN(11 , 10) , AR, LAMDA , DELTA , NX , NY , 
*  ALPHA ,CLDSRD,ELLIP,XAC,CLA 
COMMON/COF/A(101,101) ,  NEQNS 
COMMON/CAM/AA( 10 1,101) ,  AS (101 ,101) 
COMMON/FILS/IIN , IOUT , JOUT , KOUT , LOUT , MOUT , NOUT 
************************************************* 

*  * 

*  DEFINE  THE  CHARACTER  STRINGS  NEEDED  FOR  GRAPH  TITLES  * 

*  * 
*********************************************************************** 

TYPE(l)-'     FLAT' 

TYPE (2)= 'CAMBERED' 

TITLE(1)=  'd(CL)/dy 

TITLE(2)=  'd(CD)/dy 

TITLE(3)=  'd(CMAC)/dy 

TITLE(4)=  'TWIST  IN  DEGREES 

TITLE(5)=  'DELTA  CP 

TITLE (6)=  'WING  HEIGHT  COORD  Z 

TITLE2(1,1)=  ' ,  AOA  =  ' 

TITLE2(1,2)=  '  deg 

TITLE2(2,1)=  ' ,  CLi  =  ' 

TITLE2(2,2)=  ' 

AXIS(l)-   'SPAN-WISE  COORDINATE  Y 
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AXIS (2)=   'CHORD-WISE  COORDINATE  X' 
COEF(l,l)-'CL 
COEF(l,2)='CD 
COEF(l,3)-'CMAC  = 
COEF(2,l)='CLi   = 
COEF(2,2)='CD 
COEF(2,3)='CMAC  = 
**************************************************** 

*  * 

*  ASSIGN  DEFAULT  VALUES  TO  THE  VARIABLES  * 

*  * 
*********************************************************************** 

IMOD  =   1 

AR  =     2.0 

LAMDA  -1.0 

ADELTA=  1.0 

NX  =     4 

NY  =     4 

ALPHA  =5.0 

ELLIP  =  'YES' 

CLDSRD=  1.0 

ITYPE  =  1 

ISECT  =  1 

XAC  =    0.0 

ESC  -  CHAR(27) 

PI  =  2.0*ACOS(0.0) 
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****************************************************** 

*  * 

*  SET  INPUT/OUTPUT  VARIABLES  AND  OPEN  THE  NECESSARY  FILES  * 

*  * 
*********************************************************************** 

UN  =  5 

IOUT  =  6 

JOUT  =  10 

KOUT  =11 

LOUT  =12 

LDAT1  =13 

LDAT2  =  14 

OPEN(UNIT  =  UN) 

OPEN (UNIT  =  IOUT) 

OPEN (UNIT  =  JOUT,  FILE-' WING .MAT' ) 

OPEN (UNIT  =  KOUT,  FILE-' SECTION .MAT' ) 

OPEN (UNIT  =  LOUT,  FILE-' PARAMS . DAT' ) 

OPEN (UNIT  =  LDAT1,  FILE-' FLAT . DAT' ) 

OPEN (UNIT  =  LDAT2,  FILE-' CAMBER. DAT ' ) 
10    FORMAT  (I6/3(F6.2,/),2(I3/),F6.2/A/F6.2/I2/I2/F6.2) 
20    FORMAT  (/2 (10X, F12 . 6/)/10X, F12 . 6) 
*********************************************************************** 

*  * 

*  FORMAT  STATEMENTS  ARE  FOR  THE  INPUT  TO  THE  GRAPHER  PROGRAM.      * 

*  OTHER  GRAPHICS  PROGRAMS  WILL  REQUIRE  DIFFERENT  FORMAT  STMTS.     * 
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32  FORMAT  ( 

*  '1236'/ 

*  '1  2  1  0  0'/ 

*  'PI'/ 

*  '65  66  78  "NO  "  0 '  / 

*  '"NO"  "SOLID"  1.500e-001  1'/ 

*  '"YES"  41  1.000e-001  1  V / 

*  '78  9.900e+028  9.900e+028  0 . 000e+000  "DEFAULT"  1.000e-001  1'/ 

*  '"SOLID"  5  1.500e-001  0 . 000e+000  9 . 900e+029  100  2 . 000e+000  1') 

33  FORMAT  ( 

*  '1236'/ 

*  '1  2  3  0  1'/ 

*  'PI'/ 

*  '65  66  78  "NO  "  0'/ 

*  '"NO"  "SOLID"  1.500e-001  1'/ 

*  '"YES"  41  1.000e-001  1  1'/ 

*  '78  9.900e+028  9.900e+028  0.000e+000  "DEFAULT"  1.000e-001  1'/ 

*  '"SOLID"  5  1.500e-001  0 . 000e+000  9.900e+029  100  2 . 000e+000  1') 

34  FORMAT ( 

*  'HORIZ'/ 

*  '1.750e+000  ',E10.3,'  6 . 000e+000  88'/ 

*  E10.3,'  '.E10.3,'  9.900e+028  1  1'/ 

*  '0.000e+000  '.E10.3,'  1.500e-001  1  1'/ 

*  '1  1  1'/ 

*  '1  2.750e+000  4.000e-001  9.900e+028  1.800e-001'/ 
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*  ' "DEFAULT"  "DEFAULT"  " '  , A,  ""  ) 

36  FORMAT ( 

*  'Y-AXIS'/ 

*  '1.750e+000  1.000e+000  5.000e+000  89'/ 

*  E10.3,'  '.E10.3,'  9.900e+028  11'/ 

*  '2.700e+002  '.E10.3,'  1.500e-001  11'/ 

*  '1  ' ,11, '  1'/ 

*  '1  9.900e+028  0.000e+000  9.900e+028  1.800e-001'/ 

*  '"DEFAULT"  "DEFAULT"  "',A,'"') 

37  FORMAT ( 

*  'PI'/ 

*  '"DEFAULT"  1'/ 

*  '2.000e+000  7.000e+000  0 . 000e+000  1.700e-001'/ 

*  '  2 '/ 

*  '0  13  "' ,A, '  WING"'/ 

*  '1  12  "CLi  =  ' ,F5.2, '"') 

38  FORMAT ( 

*  'PI'/ 

*  '"DEFAULT"  1'/ 

*  '2.000e+000  7.000e+000  0.000e+000  1.700e-001'/ 

*  '2'/ 

*  '0  13  "' ,A, '  WING"'/ 

*  '1  34  "' ,A,F6.4,A,F5.2,A'" ') 

39  FORMAT ( 

*  'PI'/ 

*  '"DEFAULT"  1'/ 
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*  '2.000e+000  7.000e+000  O.OOOe+000  1.700e-001'/ 

*  '2'/ 

*  '0  26  "',A,'  WING,  SECTION  ',13,'"'/ 

*  '1  26  "vertical  scale  exaggerated"') 

40  FORMAT ( 

*  'PI'/ 

*  '"DEFAULT"  1'/ 

*  '2.000e+000  7.000e+000  0.000e+000  1.700e-001'/ 

*  '1'/ 

*  '0  26  " ',A,'  WING,  SECTION  ',13,'*") 

41  FORMAT ( 

*  'LE'/ 

*  '"DEFAULT"  1'/ 

*  E10.3,'  6.000e+000  0.000e+000  1.000e-001'/ 

*  '1'/ 

*  '0  4  "L.E. " ' ) 

42  FORMAT ( 

*  'TE'/ 

*  '"DEFAULT"  1'/ 

*  E10.3,'  6.000e+000  0.000e+000  1.000e-001'/ 

*  '1'/ 

*  '0  4  "T.E. "' ) 

43  FORMAT ( 

*  'PI'/ 

*  '2'/ 

*  E10.3,'  1.000e+000  ' , E10 . 3 , '  6.000e+000  1  2  1.000e-001'/ 
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*  E10.3,'  l.OOOe+000  '.E10.3,'  6 . OOOe+000  1  2  l.OOOe-001') 
50    FORMAT  (2 (F10 . 6 , IX) ) 

60  FORMAT  (23F10.6) 
70  FORMAT  (10F10.6) 
****************************************************** 

*  * 

*  THESE  FORMAT  STATEMENTS  ARE  FOR  THE  FLASHUP  WINDOWS  MENUS       * 

*  * 
*********************************************************************** 
82     FORMAT  ('  ' ,A,A,I3) 

84    FORMAT  ('  ' ,A,A,A) 
*********************************************************************** 

*  * 

*  READ  THE  LAST  SET  OF  PARAMETERS,  IF  UNCHANGED,  READ  WING  AND    * 

*  SECTN  DATA  ,  IMOD  IS  THE  FLAG  TO  TELL  IF  THE  PARAMETERS  HAVE    * 

*  BEEN  CHANGED.   1  =  CHANGE,  0  =  NO  CHANGE.  * 

*  * 

*****************************************#****#*********■*************** 

READ  (LOUT, 10)  IMOD ,AR,LAMDA, ADELTA, NX, NY, ALPHA, 

*  ELLIP,CLDSRD,ITYPE,ISECT,XAC 
DELTA  =  ADELTA*PI/180. 

NEQNS  =  NX*NY 

B  =  AR*(1.0+LAMDA)/2.0 

IF  (IMOD.EQ.O)  THEN 

READ  (JOUT.60)  ( (WING(I , J) , J=l , 23) , 1=1 , NEQNS) 

READ  (KOUT.70)  ( (SECTN(I , J) , J=l , 10) , 1=1 ,NY) 
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READ  (LDAT1.20)  (C0EFF(1 , J) , J=l , 3) 

READ  (LDAT2.20)  (C0EFF(2 , J) , J-l , 3) 
ELSE 

WRITE  (IOUT,*)  '  CHANGES  WERE  MADE  TO  THE  PLANFORM.  GO  BACK' 

WRITE  (IOUT,*)  '  AND  RUN  THE  LOAD/SHAPE  PROGRAM  AGAIN.' 

CALL  DELAY  (6) 

GO  TO  1000 
END  IF 

COEFF(l,4)  =  ALPHA 
COEFF(2,4)  =  CLDSRD 
CLOSE  (JOUT) 
CLOSE  (KOUT) 
CLOSE  (LDAT1) 
CLOSE  (LDAT2) 

*********************************************************************** 

*  * 

*  SET  MINIMUM  VALUE  OF  X.   THIS  X  IS  THE  CHORDWISE  COORDINATE 

*********************************************** 

XMIN  =0.0 
*********************************************************************** 

*  * 

*  COMPUTE  THE  CORNER   OF  THE  WING,  TO  GET  MAX  VALUES  FOR  THE  AXIS  * 

*  THIS  IS  NEEDED  IN  CASE  THE  WING  IS  SWEPT  OR  HIGHLY  TAPERED.      * 

*  * 
*********************************************************************** 
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XT  =  (B*TAN( DELTA )/2.0)  +  LAMDA 
********************************************** 

*  * 

*  NOW  CONVERT  SPAN  AND  CHORD  AXIS  LIMITS  INTO  "ROUNDED"  UNITS     * 

*  AT  THE  SAME  TIME,  COMPUTE  THE  INCREMENTS  IN  TIC  MARKS  FOR  THE    * 

*  GRAPHS .  * 

*  * 
*********************************************************************** 

IF  (XT.GT.l.)  THEN 

CALL  RNDUP(XT.XMAX) 
ELSE 

XMAX  =  1.0 
END  IF 

YMIN  =0.0 

CALL  RNDUP(B/2.0,YMAX) 
XINC  =  (XMAX-XMIN)/5.0 
YINC  =  (YMAX-YMIN)/5.0 
*********************************************************************** 

*  * 

*  CLEAR  THE  SCREEN  AND  WRITE  THE  GRAPHIC  MENU.  * 

*  * 
*********************************************************************** 
110   WRITE  (IOUT,*)  ESC,'[2J' 

WRITE  (IOUT,*)  ESC, ' [00;00H' 
WRITE  (IOUT,*)  '~C=ALL/' 
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WRITE  (IOUT,*)  '"W=GRAPH/' 

WRITE  (IOUT, 84)  ESC , ' [08 ; 69H' ,TYPE(ITYPE) 

WRITE  (IOUT, 82)  ESC , ' [09 ; 69H' , ISECT 

********************************************************************** 

*  * 

*  READ  THE  OPTION  AND  SEND  TO  THE  APPROPRIATE  AREA.  * 

*  OPTIONS  ARE:  * 

*  1,  SETS  THE  TYPE  OF  WING,  FLAT  OR  ELLIPTIC.  * 

*  2,  SETS  THE  WING  SECTION  TO  BE  PLOTTED  * 

*  3,  RETURNS  TO  THE  MAIN  MENU  * 

*  4-7  PLOT  THE  LIFT,  DRAG,  MOMENT   OR  TWIST  VS  SPAN     * 

*  8     PLOT  DELTA  CP  VERSUS  CHORD  * 

*  9      PLOTS  THE  WING  SHAPE  VERSUS  CHORD,  WHICH  GIVES     * 

*  AN  ELLIPTIC  LIFT  DISTRIBUTION.  * 

*  * 

***********************************  ************************************ 

READ  (UN,*)  ISELEC 
WRITE  (IOUT,*)  ESC, ' [ 2 J ' 
WRITE  (IOUT,*)  ESC, ' [00;00H' 
*********************************************** 

*  * 

*  THIS  IS  OPTION  1,  SET  TYPE  OF  WING.  * 

*  * 
*********************************************************************** 

IF  (ISELEC. EQ.l)  THEN 

IF  (ELLIP.EQ. 'NO  ')  THEN 
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WRITE  (IOUT,*)  '  ELLIPTIC  DATA  WAS  NOT  COMPUTED' 

CALL  DELAY(6) 

GO  TO  110 
ELSE 
ENDIF 
IF  (ITYPE.EQ.l)  THEN 

ITYPE  =  2 
ELSE 

ITYPE  =  1 
ENDIF 
GO  TO  110 
ELSE 
ENDIF 

*********************************************************************** 

*  * 

*  THIS  IS  OPTION  2,  SET  SECTION  NUMBER  * 

*************************************************** 
IF  (ISELEC.EQ.2)  THEN 

CALL  CDIT85(1,NY,ISECT, 
*   'WING  SECTION  NUMBER  ' , IMOD) 

IMOD  =  0 
GO  TO  110 
ELSE 
ENDIF 
*********************************************************************** 
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*  THIS  IS  OPTION  3,  RETURN  TO  MAIN  PROGRAM.  BEFORE  LEAVING,  WRITE  * 

*  THE  FLASHUP  WINDOW  FOR  THE  MAIN  MENU.  * 

*  * 
************************************************** 

IF  (ISELEC.EQ.3)  THEN 

GO  TO  1000 
ELSE 
ENDIF 
IF  ((ISELEC.LT.l) .OR. (ISELEC.GT.9))  THEN 

WRITE  (IOUT,*)  '  THAT  IS  NOT  A  VALID  SELECTION' 

CALL  DELAY(6) 

GO  TO  110 
ELSE 
ENDIF 
********************************************* ************************** 

*  * 

*  THESE  ARE  THE  LIFT,  DRAG  ,  MOMENT  AND  TWIST  VS  SPAN  PLOTS.       * 

*  * 
*********************************************************************** 

IF  (ISELEC.GT.7)  GO  TO  310 
*********************************************************************** 

*  * 

*  OFFSET  IS  ESTABLISHED,  IT  IS  AN  EASY  WAY  TO  USE  THE  ISELEC      * 

*  TO  FIND  THE  CORRECT  COLUMN  OF  THE  SECTN  MATRIX.  * 

*  SINCE  THERE  IS  NO  TWIST  ON  A  FLAT  WING,  SET  THE  WING  TYPE  TO     * 
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*  CAMBERED  IF  THE  TWIST  PLOT  IS  SELECTED.  * 

*  * 
*************************************************** 

IF  (ISELEC.EQ.7)  ITYPE  =  2 
IF  (ITYPE. EQ.l)  THEN 

IOFFS  =  -2 
ELSE 

IOFFS  =  3 
END  IF 
*********************************************************************** 

*  * 

*  SET  UP  THE  DATA  FILE  FOR  THE  GRAPHIC  PLOT.  * 

*  * 
*********************************************************************** 

J  =  1 

VMAX  =  -100.0 

VMIN  =  100.0 

OPEN (UNIT  =  30, FILE  ='P1.DAT') 

IJ  =  J 

DO  300  I  =  1,NY 

IJ  =  (I-1)*NX  +  J 

VERT  =  SECTN(I,IOFFS+ISELEC) 

IF  ( VMAX. LT. VERT)  VMAX  =  VERT 

IF  ( VMIN. GT. VERT)  VMIN  =  VERT 

WRITE  (30,50)  WING(IJ,5),  VERT 
300    CONTINUE 
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CLOSE  (30) 

*********************************************************************** 

*  * 

*  SET  UP  THE  Pl.GRF  FILE  FOR  GRAPHER.   IT  HAS  ALL  THE  GRAPH       * 

*  FORMATS  IN  IT.  * 

*************************************************** 
IF  (VMIN. GT. 0.0)  VMIN  =0.0 
IF  (VMAX.LT.0.0)  VMAX  =0.0 
CALL  RNDUP( VMAX, VMAX) 
CALL  RNDUP (VMIN, VMIN) 
CALL  SDIG (VMIN, VMIN, VMAX, VMAX, IDIG) 
VINC  =  (VMAX-VMIN)/5.0 
HSTART  = (VMAX- 6*VMIN)/ (VMAX- VMIN) 
OPEN (UNIT  =  32, FILE  =' Pl.GRF') 
WRITE  (32,32) 

WRITE  (32,34)  HSTART ,YMIN ,YMAX,YINC , AXIS (1) 
WRITE  (32,36)  VMIN, VMAX, VINC, IDIG, TITLE(ISELEC- 3) 
IF  (ISELEC.EQ.7)  THEN 

WRITE  (32,37)  TYPE(ITYPE) , COEFF(2 , 4) 
ELSE 

WRITE  (32,38)  TYPE(ITYPE) , COEF(ITYPE , ISELEC- 3) , 

*  COEFF(ITYPE,ISELEC-3) ,  TITLE2(ITYPE, 1) ,  COEFF(ITYPE , 4) , 

*  TITLE2(ITYPE,2) 
ENDIF 

CLOSE  (32) 

162 


GO  TO  1000 
************************************************** 

*  * 

*  THESE  ARE  THE  DELTA  CP  PLOTS  * 

*  * 
*********************************************************************** 
310    CONTINUE 

IF  (ISELEC.EQ.9)  GO  TO  400 
*********************************************************************** 

*  * 

*  OFFSET  IS  SET  AGAIN.  * 

*  * 
*********************************************************************** 

IF  (ITYPE.EQ.l)  THEN 

IOFFS  =  8 
ELSE 

IOFFS  =  10 
ENDIF 
*********************************************************************** 

*  * 

*  SET  UP  THE  PI. DAT  FILE  FOR  GRAPHER  * 

*  * 
*********************************************************************** 

I  =  ISECT 
VMAX  =  -100.0 
VMIN  =  100.0 
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OPEN (UNIT  =   30, FILE  -'PI. DAT') 
IJ   =    (I-1)*NX  +   1 
************************************************ 

*  * 

*  XLE  IS  THE  VALUE  OF  THE  LEADING  EDGE.   SINCE  THE  DELTA  CP       * 

*  FOR  THE  CAMBERED  WING  IS  ZERO  AT  THE  LEADING  AND  TRAILING       * 

*  EDGE,  IT  IS  POSSIBLE  TO  EXTEND  THE  DATA  TO  THOSE  POINTS.         * 

*  * 
*********************************************************************** 

XLE  =  WING(IJ,4)-WING(IJ,13)/2. 
IF(ITYPE.EQ.2)  THEN 

WRITE  (30,50)  XLE,  0.0 
ELSE 
END  IF 
*********************************************************************** 

*  * 

*  GET  THE  MAX  AND  MIN  VALUES  OF  THE  PARAMETER  BEING  PLOTTED       * 

*  SO  ALL  GRAPHS  WITH  THAT  PARAMETER  WILL  HAVE  THE  SAME  VERTICAL   * 

*  SCALE.  * 

*  * 
*********************************************************************** 

DO  315  IJ  =  l.NEQNS 

VERT  =  WING(IJ,IOFFS+ISELEC) 
IF  (VMAX.LT.VERT)  VMAX  =  VERT 
IF  (VMIN.GT.VERT)  VMIN  =  VERT 
315    CONTINUE 
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************************************************************ 

*  * 

*  SET  UP  THE  PI. DAT  FILE  FOR  GRAPHER  * 

*  * 
*********************************************************************** 

DO  320  J  =  1,NX 

IJ  =  (I-1)*NX  +J 
VERT  =  WING(IJ,IOFFS+ISELEC) 
WRITE  (30,50)  WING(IJ,11),  VERT 
320    CONTINUE 

IJ  =  I*NX 
*********************************************************************** 

*  * 

*  XLE  IS  COMPUTED  FOR  THE  REASONS  STATED  ABOVE. 

*  * 
*********************************************************************** 

XTE  =  WING(IJ,4)+WING(IJ,13)/2. 
IF(ITYPE.EQ.2)  THEN 

WRITE  (30,50)  XTE,  0.0 
ELSE 
END  IF 
CLOSE  (30) 
*********************************************************************** 

*  * 

*  COMPUTE  AND  STORE  THE  VALUES  FOR  THE  Pl.GRF  FILE  USED  BY  GRAPHER* 
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*************************************************** 

IF  (VMIN.GT.O.O)  VMIN  =0.0 

IF  (VMAX.LT.0.0)  VMAX  =  0.0 

CALL  RNDUP( VMAX, VMAX) 

CALL  RNDUP( VMIN, VMIN) 

CALL  SDIG (VMIN, VMIN, VMAX, VMAX, IDIG) 

VINC  =  (VMAX- VMIN) /5 . 0 

HSTART  =  (VMAX-6*VMIN)/(VMAX-VMIN) 

PLE  =  (6.*XLE/(XMAX-XMIN))+1.75 

PTE  =  (6.*XTE/(XMAX-XMIN))+1.75 

OPEN(UNIT  =  32 .FILE  ='P1.GRF') 

WRITE  (32,33) 

WRITE  (32,34)  HSTART, XMIN,XMAX,XINC , AXIS (2) 

WRITE  (32,36)  VMIN, VMAX, VINC , IDIG ,TITLE(ISELEC- 3) 

WRITE  (32,40)  TYPE(ITYPE) , ISECT 

WRITE  (32,41)  PLE 

WRITE  (32,42)  PTE 

WRITE  (32,43)  PLE , PLE , PTE , PTE 

CLOSE  (32) 

GO  TO  1000 
*********************************************************************** 

*  * 

*  THIS  IS  THE  SHAPE  PLOT.  * 

*  * 
*********************************************************************** 
400    IF  (ELLIP.EQ. 'NO  ')  THEN 
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WRITE  (IOUT,*)  '  ELLIPTIC  SHAPE  DATA  WAS  NOT  COMPUTED' 

CALL  DELAY(6) 

GO  TO  110 
ELSE 
ENDIF 
**************************************************** 

*  * 

*  SET  UP  THE  PI. DAT  FILE  FOR  THE  WING  SHAPE.   WE  KNOW  THAT  THE    * 

*  LEADING  EDGE  STARTS  AT  ZERO  HEIGHT,  SO  THAT  STARTING  VALUE  IS    * 

*  ADDED  TO  THE  DATA,  IN  ORDER  TO  GET  THE  MAXIMUM  INFORMATION      * 

*  INCLUDED  IN  THE  GRAPH  * 

*  XLE  AND  XTE  ARE  ALSO  NEEDED  SO  THAT  VERTICAL  DOTTED  LINES       * 

*  SHOWING  THEIR  POSITION  CAN  BE  ADDED  TO  THE  GRAPH.  * 

*  * 
*********************************************************************** 

I  =  ISECT 
VMAX  =  -100.0 
VMIN  -  100.0 

OPEN (UNIT  =  30, FILE  -'PI. DAT') 
IJ  =  (I-1)*NX  +  1 
XLE  =  WING(IJ,4)-WING(IJ,13)/2. 
WRITE  (30,50)  XLE,  0.0 
I J  -  I*NX 

XTE  =  WING(IJ,4)+WING(IJ,13)/2. 
DO  405  I J  =  l.NEQNS 
VERT  =  WING (IJ, 20) 
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IF  (VMAX.LT.VERT)  VMAX  =  VERT 

IF  (VMIN.GT.VERT)  VMIN  =  VERT 
405    CONTINUE 

DO  410  J  =  l.NX 

IJ  =  (I-1)*NX+J 

VERT  =  WING (I J, 20) 

WRITE  (30,50)  WING(IJ,4)+WING(IJ,13)/2.0,  VERT 
410    CONTINUE 

CLOSE  (30) 

IF  (VMIN. GT. 0.0)  VMIN  =0.0 

IF  (VMAX. LT. 0.0)  VMAX  =0.0 

CALL  RNDUP( VMAX, VMAX) 

CALL  RNDUP( VMIN, VMIN) 

CALL  SDIG (VMIN, VMIN, VMAX, VMAX, IDIG) 

VINC  =  (VMAX- VMIN) /5.0 

HSTART  =  (VMAX-6*VMIN)/(VMAX-VMIN) 

PLE  =  (6.*XLE/(XMAX-XMIN))+1.75 

PTE  =  (6.*XTE/(XMAX-XMIN))+1.75 

OPEN(UNIT  =  32, FILE  ='P1.GRF') 

WRITE  (32,33) 

WRITE  (32,34)  HSTART ,XMIN,XMAX,XINC .AXIS (2) 

WRITE  (32,36)  VMIN, VMAX, VINC , IDIG ,TITLE(ISELEC- 3) 

WRITE  (32,39)  TYPE(2) , ISECT 

WRITE  (32,41)  PLE 

WRITE  (32,42)  PTE 

WRITE  (32,43)  PLE , PLE , PTE , PTE 
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CLOSE  (32) 
1000   REWIND  (LOUT) 
*************************************************** 

*  * 

*  RE-WRITE  THE  PARAMETER  DATA  FILE,  SINCE  IT  KEEPS  TRACK  OF  THE   * 

*  LAST  WING  TYPE  AND  SECTION  PLOTTED.  * 

*  * 
*********************************************************************** 

WRITE  (LOUT, 10)  IMOD , AR , LAMDA , ADELTA , NX , NY , ALPHA , 
*ELLIP , CLDSRD , ITYPE , ISECT , XAC 
CLOSE  (LOUT) 
*********************************************************************** 

*  * 

*  BEFORE  LEAVING  THE  PROGRAM,  RE-WRITE  THE  MAIN  PROGRAM  MENU.      * 

*  * 
*********************************************************************** 

WRITE  (IOUT,*)  ESC,'[2J' 

WRITE  (IOUT,*)  ESC, ' [00;00H' 

WRITE  (IOUT,*)  '~C=ALL/' 

WRITE  (IOUT,*)  '"W-MAIN/' 

END 
$ INCLUDE: 'CDIT85.FOR' 
$ INCLUDE : ' DELAY . FOR ' 
$ INCLUDE : ' RNDUP . FOR ' 
$INCLUDE: 'SDIG.FOR' 
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DELAY 

**********************************  ************************************* 


* 
* 
* 
* 


SUBROUTINE  DELAY (INTERV)  * 

CALLED  TO  SLOW  THE  SCREEN  DOWN  AFTER  PRINTING  AN  ERROR  MESSAGE.  * 

VARIABLES  USED  ARE:  * 

INTERV  :   NUMBER  OF  SECONDS  FOR  DELAY  MUST  BE  LESS  THAN  59,  * 


* 

* 

IHR 

* 

IMIN 

* 

ISEC1 

* 

ISEC2 

* 

IHUND 

* 

INPUT. 

TIME  OF  DAY  IN  HOURS 

TIME  OF  DAY  IN  MIN 

TIME  OF  DAY  IN  SECONDS 

TIME  OF  DAY  IN  SECONDS 

TIME  OF  DAY  IN  HUNDREDTHS  OF  SECONDS. 


* 

* 
* 
* 
* 
* 


******************************************************** 

SUBROUTINE  DELAY (INTERV) 

INTEGER*2  IHR,  IMIN,  ISEC1,  ISEC2,  IHUND,  INTER 

INTER  =  INTERV 

CALL  GETTIM(IHR, IMIN, ISEC1, IHUND) 
100    CALL  GETTIM( IHR, IMIN, ISEC2, IHUND) 

IF  ((ISEC2-ISEC1) .LT. INTER)  GOTO  100 

RETURN 

END 
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CDIT85 

******************************************************* 

*  * 

*  SUBROUTINE  CDIT8 5 (LOW, HIGH, VALUE, NAME, IMOD)  * 

*  INTEGER  VERSION  OF  CDAT85  * 

*  HANDLE  INPUT  OF  PARAMETER  DATA  FOR  THE  PROGRAM. 

*  ERROR  CHECKING. 


INCLUDES  SOME  * 
* 


VARIABLES  USED  ARE: 

LOW     :   MINIMUM  VALUE  THAT  WILL  BE  ACCEPTED,  INPUT 
MAXIMUM  VALUE  THAT  WILL  BE  ACCEPTED,  INPUT 
VALUE  READ  FROM  SCREEN,  OUTPUT 
SCREEN  PROMPT  STRING,  INPUT 

FLAG  USED  IN  ANOTHER  PROGRAM,  SET  TO  1  ANYTIME 
CDAT85  IS  CALLED.,  OUTPUT 


* 

HIGH 

* 

VALUE 

* 

NAME 

* 

IMOD 

ESC 


ESCAPE  CHARACTER,  (27) 


*********************************************************************** 

SUBROUTINE  CDIT85 (LOW , HIGH , VALUE , NAME , IMOD) 

COMMON/FILS/IIN , IOUT , JOUT , KOUT , LOUT , MOUT , NOUT 

INTEGER  LOW, HIGH, VALUE 

CHARACTER* 3 5  NAME 

CHARACTERS  ESC 

ESC  =  CHAR(27) 
10    FORMAT ('  ENTER  THE  ' ,A35) 
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20  FORMAT ('  THE  MINIMUM  VALUE  FOR  ',A35,/'  IS  ',13) 

30  FORMAT ('  THE  MAXIMUM  VALUE  FOR  ',A35,/'  IS  ',13) 
*********************************************************************** 

*  * 

*  CLEAR  SCREEN,  POSITION  CURSOR,  AND  * 

*  SET  IMOD  FLAG  AND  WRITE  PROMPT  TO  SCREEN.  * 

*  * 

*********************************************************************** 

100    WRITE  (IOUT,*)  ESC,'[2J' 

WRITE  (IOUT,*)  ESC,'[10;0H' 

IMOD  =  1 

WRITE (IOUT, 10)  NAME 
*********************************************************************** 

*  * 

*  READ  SCREEN  VALUE.   ON  ERROR  OR  END  OF  DATA  SET,  GO  BACK  TO     * 

*  PROMPT .  * 

*  * 
***************************************************** 

READ(IIN,*,ERR  =  100,  END  =  100)  VALUE 
*********************************************************************** 

*  * 

*  IF  VALUES  ARE  OUT  OF  RANGE  ,  PRINT  ERROR  MSG  AND  RETURN  TO      * 

*  PROMPT .  * 

*  * 
*********************************************************************** 

IF  ( VALUE. LT. LOW)  THEN 
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WRITE  (IOUT.20)  NAME, LOW 
CALL  DELAY  (6) 
GO  TO  100 
ELSE 

IF ( VALUE. GT. HIGH)  THEN 

WRITE  (IOUT.30)  NAME, HIGH 
CALL  DELAY(6) 
GO  TO  100 
ELSE 
END  IF 
END  IF 
RETURN 
END 
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* 

HIGH 

* 

VALUE 

* 

NAME 

* 

IMOD 

CDAT85 

******************************************************** 

*  * 

*  SUBROUTINE  CDAT85 (LOW .HIGH, VALUE.NAME, IMOD)  * 

*  REAL  VERSION  * 

*  HANDLE  INPUT  OF  PARAMETER  DATA  FOR  THE  PROGRAM.  INCLUDES  SOME   * 

*  ERROR  CHECKING.  * 

*  * 

*  VARIABLES  USED  ARE:  * 

*  LOW     :   MINIMUM  VALUE  THAT  WILL  BE  ACCEPTED,  INPUT  * 

MAXIMUM  VALUE  THAT  WILL  BE  ACCEPTED,  INPUT  * 

VALUE  READ  FROM  SCREEN,  OUTPUT  * 

SCREEN  PROMPT  STRING,  INPUT  * 

FLAG  USED  IN  ANOTHER  PROGRAM,  SET  TO  1  ANYTIME  * 

*  CDAT85  IS  CALLED.,  OUTPUT  * 

*  ESC     :   ESCAPE  CHARACTER,  (27)  * 

*  * 
*********************************************************************** 

SUBROUTINE  CDAT85 (LOW , HIGH , VALUE , NAME , IMOD) 

COMMON/FILS/I IN , IOUT , JOUT , KOUT , LOUT , MOUT , NOUT 

REAL  LOW 

CHARACTER* 3 5  NAME 

CHARACTER*1  ESC 

ESC  =  CHAR(27) 
10    FORMAT ('  ENTER  THE  ' ,A35) 
20     FORMAT ('  THE  MINIMUM  VALUE  FOR  ',A35,/'  IS  \F6.2) 

174 


30    FORMAT ('  THE  MAXIMUM  VALUE  FOR  ',A35,/'  IS  ',F6.2) 
**************************************************** 

*  * 

*  CLEAR  SCREEN,  POSITION  CURSOR,  AND  * 

*  SET  IMOD  FLAG  AND  WRITE  PROMPT  TO  SCREEN.  * 

*  * 
*********************************************************************** 
100   WRITE  (IOUT,*)  ESC,'[2J' 

WRITE  (IOUT,*)  ESC,'[10;0H' 
IMOD  =  1 

WRITE (IOUT, 10)  NAME 
*********************************************************************** 

*  * 

*  READ  SCREEN  VALUE.   ON  ERROR  OR  END  OF  DATA  SET,  GO  BACK  TO     * 

*  PROMPT .  * 

*  * 
*********************************************************************** 

READ  (UN,*,  ERR  =  100,  END  =  100)  VALUE 
*********************************************************************** 

*  * 

*  IF  VALUES  ARE  OUT  OF  RANGE  ,  PRINT  ERROR  MSG  AND  RETURN  TO      * 

*  PROMPT .  * 

*  * 
*********************************************************************** 


IF  ( VALUE. LT. LOW)  THEN 

WRITE  (IOUT, 20)  NAME , LOW 
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CALL  DELAY  (6) 
GO  TO  100 
ELSE 

IF(VALUE.GT.HIGH)  THEN 

WRITE  (IOUT.30)  NAME .HIGH 
CALL  DELAY  (6) 
GO  TO  100 
ELSE 
END  IF 
END  IF 
RETURN 
END 


176 


GAUSS 

*********************************************************************** 

*  * 

*  SUBROUTINE  GAUSS (NRHS)  * 

*  * 

*  SOLUTION  OF  LINEAR  ALGEGRAIC  SYSTEM  BY  GAUSS  ELIMINATION  * 

*  WITH  PARTIAL  PIVOTING.   TAKEN  FROM  'AN  INTRODUCTION  TO  * 

*  THEORETICAL  AND  COMPUTATIONAL  AERODYNAMICS',  JACK  MORAN,  * 

*  JOHN  WILEY  &  SONS,  NEW  YORK,  1984,  PG  78.  * 

*  * 

*  VARIABLES  USED  ARE:  * 

*  A()     :   COEFFICIENT  MATRIX,  AUGMENTED  WITH  RIGHT  HAND  * 

*  SIDES  IN  COLUMN  NEQNS+1 ,  INPUT/OUTPUT,  UNUSABLE  * 

*  AFTER  RETURN  * 

*  NEQNS   :   NUMBER  OF  EQUATIONS  IN  MATRIX,  INPUT  * 

*  NRHS    :   NUMBER  OF  RIGHT  HAND  SIDES,  INPUT  * 

*  * 
*********************************************************************** 

SUBROUTINE  GAUSS (NRHS) 
COMMON/COF/A(101,101) , NEQNS 
NP  =  NEQNS+1 
NTOT  =  NEQNS+NRHS 
*********************************************************************** 

*  * 

*  SEARCH  FOR  THE  LARGEST  ENTRY  IN  THE  (I-l)TH  COLUMN,  ON  OR  * 

*  BELOW  THE  MAIN  DIAGONAL.  * 
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*********************************************************************** 

DO  150  I  =  2.NEQNS 
IM  =  1-1 
IMAX  =  IM 

AMAX  =  ABS(A(IM,IM)) 
DO  110  J  =  I.NEQNS 

IF  (AMAX.GE.ABS(A(J,IM)))  GO  TO  110 
IMAX  =  J 

AMAX  =  ABS(A(J,IM)) 
110       CONTINUE 
******************************************************* 

*  * 

*  SWITCH  THE  (I-l)TH  AND  IMAXTH  EQUATIONS  * 

*  * 
*********************************************************************** 

IF  (IMAX.NE.IM)  GO  TO  140 
DO  130  J  =  IM.NTOT 
TEMP  =  A(IM,J) 
A(IM,J)  =  A(IMAX,J) 
A(IMAX,J)  =  TEMP 
130      CONTINUE 
*********************************************************************** 

*  * 

*  ELIMINATE  THE  (I-l)TH  UNKNOWN  FROM  ITH  THROUGH  (NEQNS)TH        * 

*  EQUATIONS . 
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********************************************************* 
140      DO  150  J  =  I, NEONS 

R  =  A(J,IM)/A(IM,IM) 

DO  150  K  =  I.NTOT 

A(J,K)  -  A(J,K)-R*A(IM,K) 
150    CONTINUE 
*********************************************************************** 

*  * 

*  BACK  SUBSTITUTE  * 

*  * 
*********************************************************************** 

DO  220  K  =  NP,NTOT 

A(NEQNS.K)  =  A(NEQNS,K)/A(NEQNS,NEQNS) 
DO  210  L  =  2.NEQNS 
I  =  NEQNS+1-L 
IP  =  1+1 
DO  200  J  =  IP.NEQNS 

A(I,K)  =  A(I,K)  -  A(I,J)*A(J,K) 
200  CONTINUE 

A(I,K)  =  A(I,K)/A(I,I) 
210      CONTINUE 
220    CONTINUE 
RETURN 
END 
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RNDUP 

*********************************************************************** 

*  * 

*  SUBROUTINE  TO  ROUND  UP  THE  VALUE  OF  AN  AXIS  LIMIT.  * 

*  FOR  EXAMPLE  IF  THE  VALUE  PASSED  IS  .0018,  THE  VALUE  RETURNED  * 

*  IS  .0020.   THIS  IS  NECESSARY  TO  ENSURE  THAT  AXIS  LIMITS  ARE  * 

*  WELL  DEFINED.  * 

*  * 

*  VARIABLES  USED  ARE:  * 

*  VALIN    THE  VALUE  PASSED  * 

*  VALOUT   THE  VALUE  RETURNED  * 

*  VAL     A  TEMPORARY  VARIABLE  * 

*  IRTN     AN  INTEGER  MULTIPLIER  USED  TO  RETURN  TO  THE  CORRECT  * 

*  ORDER  OF  MAGNITUDE.  * 

*  ISIDE    AN  INTEGER  ADDED  TO  ROUND  UP.   DEPENDS  ON  WHETHER  THE    * 

*  VALUE  PASSED  IS  POSITIVE  OR  NEGATIVE.  * 

*  IV      A  TEMPORARY  VARIABLE  * 

******************************************************** 
SUBROUTINE  RNDUP (VALIN, VALOUT) 
RTN  =1.0 
VAL  =  VALIN 
IF  (VAL. EQ. 0.0)  THEN 
VALOUT  =  VAL 
GOTO  200 
ELSE 
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END  IF 

IF  (VAL.LT.O.O)  THEN 

ISIDE  =  -1 
ELSE 

ISIDE  =  1 
ENDIF 
100    IF  (ABS(VAL) .GE.10.0)  THEN 
VAL  =  VAL/10.0 
RTN  =  RTN*10.0 
ELSE 

IF(ABS(VAL).LT.1.0)  THEN 
VAL  =  VAL*10.0 
RTN  =  RTN/10.0 
ELSE 

IV  =  AINT(VAL)+ISIDE 
VALOUT  =  REAL(IV*RTN) 
GO  TO  200 
ENDIF 
ENDIF 
GO  TO  100 
200   RETURN 
END 
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SDIG 

*************************************************** 

*  * 

*  SUBROUTINE  SDIG(MININ,MINOUT , MAXIN , MAXOUT , IDIG)  * 

*  * 

*  SUBROUTINE  TO  GET  THE  RANGE  OF  VALUES  IN  THE  AXIS  LIMITS        * 

*  TO  BE  5  UNITS  APART.  AND  RETURN  THE  NUMBER  OF  SIGNIFICANT  * 

*  DIGITS  FOR  THE  AXIS.  * 

*  * 

*  VARIABLES  USED  ARE:  * 

*  MININ    THE  MINIMUM  AXIS  VALUE  PASSED  * 

*  MINOUT  THE  MINIMUM  AXIS  VALUE  RETURNED  * 

*  MAXIN    THE  MAXIMUM  AXIS  VALUE  PASSED  * 

*  MAXOUT   THE  MAXIMUM  AXIS  VALUE  RETURNED  * 

*  VAL     A  TEMPORARY  VARIABLE  * 

*  VALNOT   A  TEMPORARY  VARIABLE  * 

*  RTN     A  MULTIPLIER  USED  TO  RETURN  TO  THE  CORRECT  * 

*  ORDER  OF  MAGNITUDE.  * 

*  ISIDE   AN  INTEGER  TO  DETERMINE  WHETHER  THE  VALUE  IS  POSITIVE   * 

*  OR  NEGATIVE  * 

*  IV      A  TEMPORARY  VARIABLE  * 

*  * 
*********************************************************************** 

SUBROUTINE  SDIG (MININ , MINOUT , MAXIN , MAXOUT , IDIG) 
REAL  MININ , MINOUT , MAXIN , MAXOUT , RTN , VAL , VALNOT 
INTEGER  IDIG , ICASE , ISIDE , IV , ITEST 
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RTN   =1.0 

IF(ABS(MININ) .LT.ABS(MAXIN))  THEN 
************************************************** 

*  * 

*  ICASE  -  1  IS  FOR  THE  MIN  AXIS  VALUE  BEING  SMALLER  IN  MAGNITUDE   * 

*  THAN  THE  MAX  AXIS  VALUE.  * 

*  * 
*********************************************************************** 

ICASE  =  1 
VAL  =  MAXIN 
VALNOT  =  MININ 
ELSE 
*********************************************************************** 

*  * 

*  ICASE  =  2  IS  FOR  THE  MAX  AXIS  VALUE  BEING  SMALLER  IN  MAGNITUDE   * 

*  THAN  THE  MIN  AXIS  VALUE.  * 

*  * 
*********************************************************************** 

ICASE  =  2 
VAL  =  MININ 
VALNOT  =  MAXIN 
ENDIF 
*********************************************************************** 

*  * 

*  DETERMINE  WHETHER  THE  SMALLER  (IN  MAGNITUDE)  AXIS  VALUE  IS      * 

*  POSITIVE  OR  NEGATIVE.  * 
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**************************************************** 

ISIDE  =  SIGN(l.VALNOT) 
100    IF  (ABS(VAL) .GE.10.0)  THEN 
VAL  =  VAL/10.0 
RTN  =  RTN*10.0 
ELSE 

IF(ABS(VAL) .LT.1.0)  THEN 
VAL  =  VAL*10.0 
RTN  =  RTN/10.0 
ELSE 

IF  (VALNOT.EQ.0.0)  GOTO  300 
IV  =  NINT(VALNOT/RTN) 
IF(IV.EQ.O)  IV  =  IV  +  ISIDE 
GO  TO  200 
END  IF 
ENDIF 
GO  TO  100 
200    CONTINUE 

ITEST  =  ABS(NINT(VAL)-IV) 

IF  ((ITEST. EQ. 5) .OR. (ITEST . EQ. 10) .OR. (ITEST. EQ. 15) .OR. 

*  (ITEST. EQ. 20) .OR. (NINT(VAL) .EQ.-IV))  GOTO  300 
IV  -  IV  +  ISIDE 

ITEST  =  ABS(NINT(VAL)-IV) 

IF  ((ITEST. EQ. 5) .OR. (ITEST. EQ. 10) .OR. (ITEST . EQ. 15) .OR. 

*  (ITEST. EQ. 20) .OR. (NINT(VAL) .EQ. -IV))  GOTO  300 
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VAL  =  VAL  -  I SIDE 

ITEST  =  ABS(NINT(VAL)-IV) 

IF  ((ITEST. EQ. 5) .OR. (ITEST. EQ. 10) .OR. (ITEST . EQ. 15) .OR. 

*  (ITEST. EQ. 20) .OR. (NINT(VAL) .EQ.-IV))  GOTO  300 
IV  =  IV  +  ISIDE 

ITEST  =  ABS(NINT(VAL)-IV) 

IF  ((ITEST. EQ. 5) .OR. (ITEST. EQ. 10) .OR. (ITEST . EQ. 15) .OR. 

*  ( ITEST. EQ. 20) .OR. (NINT(VAL) .EQ.-IV))  GOTO  300 
VAL  =  VAL  -  ISIDE 

300    IDIG  =  NINT(LOG10(1.0/RTN))  +  1 
IF  (ICASE.EQ.l)  THEN 

MINOUT  =  IV*RTN 

MAXOUT  =  NINT(VAL)*RTN 
ELSE 

MAXOUT  =  IV*RTN 

MINOUT  =  NINT(VAL)*RTN 
END  IF 
RETURN 
END 
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APPENDIX  B.  NOMENCLATURE 

AR  aspect  ratio 

b  span 

c  local  chord 

c  average  chord 

c  r  root  chord 

c  t  tip  chord 

Cj)  drag  coefficient 

Cj)i  induced  drag  coefficient 

C £  lift  coefficient 

Cj^i  ideal  lift  coefficient 

Cm  moment   coefficient   about   the   leading   edge   of   an 

airfoil 

Cm  root  chord  for  circulation  model 

Cfl  moment  coefficient 

Ctf0  moment  coefficient  about  the  y  axis 

^Mac  moment  coefficient  about  the  aerodynamic  center 

ACp  pressure  difference  coefficient 

d  grid  element  inset 

D  drag 

AD i  section  induced  drag 

AF  force  per  unit  span 

K  vortex  drag  factor 

L  lift 
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s 
vYff 

w 


x 


ac 


cp 


M  moment 

r         distance  between  field  and  control  points 

planform  area 

local  effective  velocity- 
remote  velocity 

induced  velocity 

chordwise  coordinate 

x  coordinate  of  aerodynamic  center 

chordwise  position  of  center  of  pressure 

spanwise  coordinate 

height  coordinate 

field  point 
(xO'Yo)     control  point 

angle  of  attack 

c  ircula t ion 

incremental  circulation 

leading  edge  sweep  angle 

transformed  spanwise  coordinate 

tangent  of  leading  edge  sweep  angle  in  circulation 

mode  1 
A        taper  ratio  in  VORTEX  model 
p        density 

r        percent  wing  taper  in  circulation  model 
<f>  transformed  chordwise  coordinate 
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a 
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control  po  int 
field  point 
wing  shape 


position  where  velocity  is  evaluated 
position  of  the  vortex  element 
surface  slope  of  the  wing,   including  twist 
from  root  to  tip  and  camber 
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